Solving high-dimensional optimal stopping problems using deep learning

TL;DR

This paper introduces a deep learning-based algorithm to efficiently solve high-dimensional optimal stopping problems, particularly for pricing complex American and Bermudan options with up to 5000 dimensions, overcoming the curse of dimensionality.

Contribution

The paper presents a novel deep learning algorithm that approximates optimal exercise strategies and prices for high-dimensional optimal stopping problems, applicable beyond finance.

Findings

Effective in high dimensions up to 5000 assets

Achieves high accuracy compared to reference values

Demonstrates fast computation times

Abstract

Nowadays many financial derivatives, such as American or Bermudan options, are of early exercise type. Often the pricing of early exercise options gives rise to high-dimensional optimal stopping problems, since the dimension corresponds to the number of underlying assets. High-dimensional optimal stopping problems are, however, notoriously difficult to solve due to the well-known curse of dimensionality. In this work, we propose an algorithm for solving such problems, which is based on deep learning and computes, in the context of early exercise option pricing, both approximations of an optimal exercise strategy and the price of the considered option. The proposed algorithm can also be applied to optimal stopping problems that arise in other areas where the underlying stochastic process can be efficiently simulated. We present numerical results for a large number of example problems, which include the pricing of many high-dimensional American and Bermudan options, such as Bermudan max-call options in up to 5000 dimensions. Most of the obtained results are compared to reference values computed by exploiting the specific problem design or, where available, to reference values from the literature. These numerical results suggest that the proposed algorithm is highly effective in the case of many underlyings, in terms of both accuracy and speed.

Tables

Table 1. Table 1: Numerical simulations of the algorithm of Framework  3.2 for optimally stopping a correlated Brownian motion in the case of the Bermudan put option with three exercise opportunities of Subsection  4.3.1.1 with H = 2 𝐻 2 H=2 . In the approximative calculations of the relative approximation errors, the exact number ( 100 ) was replaced by the value 7.894 7.894 7.894 , which was approximatively computed in Matlab .

Dimen-	Mean	Standard	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	approx.	deviation	runtime
		of $𝒫$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $𝒫$
1	7.896	0.005	0.0005	0.0003	2.9
5	7.893	0.007	0.0007	0.0005	2.8
10	7.895	0.004	0.0005	0.0003	2.7
50	7.889	0.005	0.0007	0.0005	2.8
100	7.890	0.002	0.0005	0.0002	3.0
500	7.894	0.005	0.0004	0.0003	7.8
1000	7.892	0.005	0.0006	0.0004	17.2

Table 2. Table 2: Numerical approximations of the Bermudan put option with three exercise opportunities of Subsection  4.3.1.1 based on logistic regressions using only the pay-off as input ( H = 0 𝐻 0 H=0 ).

Dimen-	Mean	Standard	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	approx.	deviation	runtime
		of $𝒫$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $𝒫$
1	7.893	0.007	0.0006	0.0007	1.6
5	7.896	0.005	0.0005	0.0003	1.6
10	7.895	0.005	0.0005	0.0003	1.6
50	7.892	0.004	0.0005	0.0003	1.7
100	7.890	0.004	0.0005	0.0004	1.8
500	7.893	0.005	0.0005	0.0005	3.0
1000	7.894	0.003	0.0003	0.0002	5.7

Table 3. Table 3: Approximations of the price of the three-exercise-opportunities put option of Subsection  4.3.1.1 for d = 10 𝑑 10 d=10 with networks with ( d + 1 ) 𝑑 1 (d+1) -dimensional input layers and H ∈ { 0 , 1 , 2 , 3 , 4 , 5 } 𝐻 0 1 2 3 4 5 H\in\{0,1,2,3,4,5\} hidden layers.

Number	Mean	Standard	Rel. $L^1$-	Standard	Average
of hidden	of $𝒫$	deviation	approx.	deviation	runtime
layers $H$		of $𝒫$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $𝒫$
0	7.894	0.005	0.0005	0.0002	1.9
1	7.892	0.005	0.0006	0.0004	1.9
2	7.892	0.006	0.0006	0.0004	1.8
3	7.891	0.004	0.0005	0.0003	1.7
4	7.888	0.009	0.0007	0.0011	1.7
5	7.889	0.006	0.0008	0.0007	1.8

Table 4. Table 4: Numerical simulations of the algorithm of Framework  3.2 for optimally stopping a standard Brownian motion in the case of the American put-type example in Subsection  4.3.1.2 . In the approximative calculations of the relative approximation errors, the exact number ( 103 ) was replaced by the value 5.318 5.318 5.318 , which was obtained using the binomial tree method on Smirnov’s website  [ 87 ] .

Dimen-	Mean	Standard	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	approx.	deviation	runtime
		of $𝒫$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $𝒫$
1	5.311	0.004	0.0014	0.0007	37.9
5	5.310	0.003	0.0015	0.0006	43.2
10	5.310	0.003	0.0016	0.0006	44.4
50	5.307	0.003	0.0020	0.0005	75.0
100	5.305	0.004	0.0024	0.0007	73.8
500	5.299	0.004	0.0035	0.0007	347.7
1000	5.296	0.005	0.0042	0.0009	696.7

Table 5. Table 5: Numerical simulations of the algorithm of Framework  3.2 for pricing the American geometric average put-type option from the example in Subsection  4.3.2.1 . In the approximative calculations of the relative approximation errors, the exact value of the price ( 108 ) was replaced by the value 6.545 6.545 6.545 , which was obtained using the binomial tree method on Smirnov’s website  [ 87 ] .

Dimen-	Mean	Standard	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	approx.	deviation	runtime
		of $𝒫$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $𝒫$
40	6.512	0.004	0.0051	0.0006	166.4
80	6.509	0.003	0.0056	0.0005	263.8
120	6.507	0.003	0.0058	0.0004	350.9
160	6.504	0.003	0.0062	0.0005	358.5
200	6.501	0.006	0.0067	0.0009	432.2

Table 6. Table 6: Numerical simulations of the algorithm from  [ 9 ] for pricing the American geometric average put-type option from the example in Subsection  4.3.2.1 . In the approximative calculations of the relative approximation errors, the exact value of the price ( 108 ) was again replaced by the value 6.545 6.545 6.545 .

Dimen-	Mean	Standard	Rel. $L^1$-	Standard	Average
sion $d$	of $\hat{L}$	deviation	approx.	deviation	runtime
		of $\hat{L}$	error	of the rel.	in sec. for
				approx.	one realisa-
				error	tion of $\hat{L}$
40	6.505	0.003	0.0061	0.0005	167.8
80	6.503	0.003	0.0065	0.0005	264.9
120	6.502	0.003	0.0066	0.0004	381.1
160	6.505	0.006	0.0061	0.0008	359.7
200	6.506	0.005	0.0060	0.0008	433.2

Table 7. Table 7: Numerical simulations of the algorithm of Framework  3.2 for pricing the American geometric average call-type option from the example in Subsection  4.3.2.2 . In the approximative calculations of the relative approximation errors, the exact value of the price ( 113 ) was replaced by the number ( 114 ), which was approximatively calculated using the binomial tree method on Smirnov’s website  [ 87 ] .

Dimen-	Mean	Standard	Price	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	(114)	approx.	deviation	runtime
		of $𝒫$		error	of the rel.	in sec. for
					approx.	one realisa-
					error	tion of $𝒫$
3	0.10698	0.00006	0.10719	0.0020	0.0005	41.7
20	0.10006	0.00006	0.10033	0.0027	0.0006	57.0
100	0.09905	0.00006	0.09935	0.0030	0.0006	137.2
200	0.09891	0.00009	0.09923	0.0032	0.0009	243.0

Table 8. Table 8: Numerical simulations of the algorithm of Framework  3.2 for pricing the American geometric average call-type option from the example in Subsection  4.3.2.3 . In the approximative calculations of the relative approximation errors, the exact value of the price ( 118 ) was replaced by the number ( 119 ), which was approximatively computed in Matlab .

Dimen-	Mean	Standard	Price	Rel. $L^1$-	Standard	Average
sion $d$	of $𝒫$	deviation	(119)	approx.	deviation	runtime
		of $𝒫$		error	of the rel.	in sec. for
					approx.	one realisa-
					error	tion of $𝒫$
40	23.6878	0.0037	23.6883	0.00012	0.00009	70.8
80	23.7221	0.0019	23.7235	0.00008	0.00006	112.5
120	23.7361	0.0028	23.7357	0.00011	0.00004	147.9
160	23.7414	0.0017	23.7419	0.00006	0.00004	188.3
200	23.7454	0.0017	23.7456	0.00005	0.00004	224.5
400	23.7523	0.0010	23.7531	0.00004	0.00003	1088.3

Table 9. Table 9: Numerical simulations of the algorithm of Framework  3.2 for pricing the Bermudan max-call option from the example in Subsection  4.4.1.1 for d ∈ { 2 , 3 , 5 } 𝑑 2 3 5 d\in\{2,3,5\} and H = 2 𝐻 2 H=2 (cf.  Python code  2 in Subsection  5.2 below).

Di-	Ini-	Mean	Standard	Bino-	95% confidence	95% confidence	Average
men-	tial	of $𝒫$	deviation	mial	interval in [3]	interval in [24]	runtime
sion	value		of $𝒫$	value			in sec. for
$d$	${\xi }_1$			in [3]			one realisa-
							tion of $𝒫$
2	90	8.068	0.006	8.075	[8.053, 8.082]	–	20.9
2	100	13.901	0.010	13.902	[13.892, 13.934]	–	21.0
2	110	21.341	0.008	21.345	[21.316, 21.359]	–	21.0
3	90	11.278	0.010	11.29	[11.265, 11.308]	–	21.2
3	100	18.685	0.006	18.69	[18.661, 18.728]	–	21.3
3	110	27.560	0.007	27.58	[27.512, 27.663]	–	21.3
5	90	16.631	0.009	–	[16.602, 16.655]	[16.620, 16.653]	21.6
5	100	26.147	0.009	–	[26.109, 26.292]	[26.115, 26.164]	21.7
5	110	36.774	0.009	–	[36.704, 36.832]	[36.710, 36.798]	21.7

Table 10. Table 10: Numerical simulations of the algorithm of Framework  3.2 for pricing the Bermudan max-call option from the example in Subsection  4.4.1.1 for d ∈ { 10 , 20 , 30 , 50 , 100 , 200 , 500 } 𝑑 10 20 30 50 100 200 500 d\in\{10,20,30,50,\allowbreak 100,200,500\} and H = 2 𝐻 2 H=2 (cf.  Python code  2 in Subsection  5.2 below).

Dimen-	Initial	Mean	Standard	Average runtime
sion $d$	value ${\xi }_1$	of $𝒫$	deviation	in sec. for one
			of $𝒫$	realisation of $𝒫$
10	90	26.196	0.012	22.7
10	100	38.272	0.008	22.9
10	110	50.812	0.007	22.9
20	90	37.692	0.011	26.9
20	100	51.572	0.010	26.9
20	110	65.510	0.010	27.5
30	90	44.831	0.014	30.3
30	100	59.516	0.010	30.0
30	110	74.234	0.009	30.0
50	90	53.897	0.016	36.8
50	100	69.572	0.008	36.8
50	110	85.262	0.011	36.8
100	90	66.359	0.010	55.6
100	100	83.390	0.014	55.6
100	110	100.421	0.014	55.6
200	90	79.005	0.009	94.3
200	100	97.414	0.012	94.4
200	110	115.832	0.012	94.4
500	90	95.970	0.012	265.0
500	100	116.254	0.010	261.4
500	110	136.534	0.015	264.9

Table 11. Table 11: Numerical results for the Bermudan max-call option from the example in Subsection  4.4.1.1 for d = 5 𝑑 5 d=5 and ξ 1 = 100 subscript 𝜉 1 100 \xi_{1}=100 obtained using networks with H ∈ { 0 , 1 , 2 , 3 , 4 , 5 } 𝐻 0 1 2 3 4 5 H\in\{0,1,2,3,4,5\} hidden layers.

Number	Mean	Standard	Average runtime
of hidden	of $𝒫$	deviation	in sec. for one
layers $H$		of $𝒫$	realisation of $𝒫$
0	25.843	0.010	11.4
1	26.134	0.009	21.7
2	26.147	0.007	37.7
3	26.147	0.012	54.6
4	26.146	0.006	70.2
5	26.144	0.008	84.7

Table 12. Table 12: Numerical simulations of the algorithm of Framework  3.2 for pricing the Bermudan max-call option on 5000 stocks from the example in Subsection  4.4.1.2 (cf.  Python code  2 in Subsection  5.2 below). In the approximative calculations of the relative approximation error, the exact value of the price ( 126 ) was replaced by the value 165.430, which corresponds to a realisation of 𝒫 𝒫 \mathcal{P} with M = 6000 𝑀 6000 M=6000 .

Number of	Realisation	95% confidence	Rel. approx.	Runtime
steps $M$	of $𝒫$	interval	error	in sec.
0	106.711	[106.681, 106.741]	0.35495	157.3
250	132.261	[132.170, 132.353]	0.20050	271.7
500	156.038	[155.975, 156.101]	0.05677	386.0
750	103.764	[103.648, 103.879]	0.37276	500.4
1000	161.128	[161.065, 161.191]	0.02601	614.3
1250	162.756	[162.696, 162.816]	0.01616	728.8
1500	164.498	[164.444, 164.552]	0.00563	842.8
1750	163.858	[163.803, 163.913]	0.00950	957.3
2000	165.452	[165.400, 165.505]	0.00014	1071.9
6000	165.430	[165.378, 165.483]	0.00000	2899.5

Table 13. Table 13: Numerical simulations of the algorithm of Framework  3.2 for pricing the Bermudan max-call option from the example in Subsection  4.4.1.3 (cf.  Python code  3 in Subsection  5.3 below).

Dimen-	Matu-	Strike	Mean	Standard	European	Price	Average runtime
sion $d$	rity $T$	price $K$	of $𝒫$	deviation	price (131)	in [5]	in sec. for one
				of $𝒫$			realisation of $𝒫$
10	1	35	10.364	0.001	10.365	10.36	15.8
10	1	40	5.540	0.002	5.540	5.54	15.8
10	1	45	1.895	0.001	1.896	1.90	15.8
10	4	35	16.518	0.004	16.520	16.53	15.7
10	4	40	11.867	0.004	11.870	11.87	15.6
10	4	45	7.799	0.005	7.804	7.81	15.6
10	7	35	20.914	0.006	20.916	20.92	15.6
10	7	40	16.373	0.009	16.374	16.38	15.5
10	7	45	12.270	0.006	12.277	12.28	15.5
400	12	35	55.714	0.012	55.714	–	126.7
400	12	40	50.967	0.012	50.964	–	127.6
400	12	45	46.235	0.012	46.234	–	126.4

Table 14. Table 14: Numerical simulations of the algorithm of Framework  3.2 for pricing the American strangle spread basket option from the example in Subsection  4.4.2 (cf.  Python code  4 in Subsection  5.4 below).

Mean	Standard	Lower	Average runtime
of $𝒫$	deviation	bound	in sec. for one
	of $𝒫$	in [66]	realisation of $𝒫$
11.797	0.004	11.75	21.7

Table 15. Table 15: Numerical simulations of the algorithm of Framework  3.2 for pricing the American put basket option in Dupire’s local volatility model from the example in Subsection  4.4.3 (cf.  Python code  5 in Subsection  5.5 below). The corresponding European put basket option price is approximately equal to the value 1.741 1.741 1.741 in the case δ = 0 % 𝛿 percent 0 \delta=0\% and to the value 6.304 6.304 6.304 in the case δ = 10 % 𝛿 percent 10 \delta=10\% .

Divi-	Time discreti-	Mean	Standard	Average runtime
dend	sation para-	of $𝒫$	deviation	in sec. for one
yield $\delta$	meter $N$		of $𝒫$	realisation of $𝒫$
0%	5	1.934	0.001	8.7
0%	10	1.977	0.001	12.1
0%	50	1.976	0.001	34.2
0%	100	1.972	0.002	67.0
10%	5	6.301	0.003	8.2
10%	10	6.303	0.004	12.0
10%	50	6.305	0.002	34.2
10%	100	6.303	0.003	67.0

Table 16. Table 16: Numerical simulations of the algorithm of Framework  3.2 for pricing the path-dependent financial derivative from the example in Subsection  4.4.4 (cf.  Python code  6 in Subsection  5.6 below). According to  [ 90 , Subsection IV.D] , the value 1.282 1.282 1.282 is a lower bound for the price ( 145 ).

Time	Mean	Standard	Average runtime
horizon $T$	of $𝒫$	deviation	in sec. for one
		of $𝒫$	realisation of $𝒫$
100	1.2721	0.0001	212.2
150	1.2820	0.0002	683.7
200	1.2895	0.0002	293.0
250	1.2959	0.0001	408.4
1000	1.2998	0.0007	864.5

Equations

$$d{X}_t=\mu (X_t)dt+\sigma (X_t)d{W}_t,X_0=\xi ,t\in [0,T],$$

$$\sup\!\left\{\mathbb{E}\big{[}g(\tau,X_{\tau})\big{]}\colon\begin{subarray}{c}\tau\colon\Omega\to[0,T]\text{ is an}\\ \mathbb{F}\text{-stopping time}\end{subarray}\right\}.$$

$$0=t_0<t_1<\cdots <t_N=T$$

$$\begin{array}{rl}{X}_{t_{n+1}}& =X_{t_n}+{\int }_{t_n}^{t_{n+1}}\mu (X_s)ds+{\int }_{t_n}^{t_{n+1}}\sigma (X_s)d{W}_s.\end{array}$$

$$\begin{array}{rl}{X}_{t_{n+1}}& \approx X_{t_n}+\mu (X_{t_n})\left(t_{n+1}-t_n\right)+\sigma (X_{t_n})\left(W_{t_{n+1}}-W_{t_n}\right).\end{array}$$

$${\mathcal{X}}_{n+1}={\mathcal{X}}_n+\mu ({\mathcal{X}}_n)\left(t_{n+1}-t_n\right)+\sigma ({\mathcal{X}}_n)\left(W_{t_{n+1}}-W_{t_n}\right)$$

$$\begin{split}\sup\!\left\{\mathbb{E}\big{[}g(t_{\tau},\mathcal{X}_{\tau})\big{]}\colon\begin{subarray}{c}\tau\colon\Omega\to\{0,1,\dots,N\}\text{ is an}\\ \mathfrak{F}\text{-stopping time}\end{subarray}\right\}\approx\sup\!\left\{\mathbb{E}\big{[}g(\tau,X_{\tau})\big{]}\colon\begin{subarray}{c}\tau\colon\Omega\to[0,T]\text{ is an}\\ \mathbb{F}\text{-stopping time}\end{subarray}\right\}\end{split}$$

$$Y=f\circ X.$$

$$\Omega\ni\omega\mapsto\sum_{n=0}^{N}n\,\mathbb{U}_{n}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{n}(\omega)\big{)}\in\{0,1,\dots,N\}$$

$$\tau =\sum _{n=0}^{N}n{\mathbb{U}}_n({\mathcal{Z}}_0,{\mathcal{Z}}_1,\dots ,{\mathcal{Z}}_n).$$

$$\begin{split}&\left\{\omega\in\Omega\colon\sum_{n=0}^{N}n\,\mathbb{U}_{n}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{n}(\omega)\big{)}=k\right\}\\ &=\left\{\omega\in\Omega\colon\mathbb{U}_{k}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{k}(\omega)\big{)}=1\right\}\\ &=\left\{\omega\in\Omega\colon(\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{k}(\omega))\in\underbrace{(\mathbb{U}_{k})^{-1}(\{1\})}_{\in\mathcal{B}((\mathbb{R}^{d+1})^{k+1})}\right\}\in\mathfrak{F}_{k}.\end{split}$$

$$\varrho (\omega )=\sum _{n=0}^{N}n\text{mathbbm}{1}_{\{\varrho =n\}}(\omega ).$$

$$\Omega \ni \omega \mapsto \text{mathbbm}{1}_{\{\tau =n\}}(\omega )\in \{0,1\}$$

$$\forall\,n\in\{0,1,\dots,N\}\colon\sigma_{\Omega}\big{(}(\mathcal{Z}_{0},\mathcal{Z}_{1},\dots,\mathcal{Z}_{n})\big{)}=\mathfrak{F}_{n}$$

$$\Omega \ni \omega \mapsto \text{mathbbm}{1}_{\{\tau =n\}}(\omega )\in \{0,1\}$$

$$\mathbbm{1}_{\{\tau=n\}}(\omega)=\mathbb{V}_{n}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{n}(\omega)\big{)}.$$

$$\begin{array}{rl}& {\mathbb{U}}_n(z_0,z_1,\dots ,z_n)\\ & =\max \left\{{\mathbb{V}}_n(z_0,z_1,\dots ,z_n),n+1-N\right\}\left[1-\sum _{k=0}^{n-1}{\mathbb{U}}_k(z_0,z_1,\dots ,z_k)\right].\end{array}$$

$${\mathbb{U}}_N(z_0,z_1,\dots ,z_N)=\left[1-\sum _{k=0}^{N-1}{\mathbb{U}}_k(z_0,z_1,\dots ,z_k)\right].$$

$$\sum _{k=0}^N{\mathbb{U}}_k(z_0,z_1,\dots ,z_k)=1.$$

$${\mathbb{U}}_0(z_0)={\mathbb{V}}_0(z_0).$$

$$\forall n\in \{0,1,\dots ,N\},z_0,z_1,\dots ,z_n\in {\mathbb{R}}^{d+1}:{\mathbb{V}}_n(z_0,z_1,\dots ,z_n)\in \{0,1\},$$

$$\left\{{\mathbb{U}}_0(z_0),{\mathbb{U}}_1(z_0,z_1),\dots ,{\mathbb{U}}_n(z_0,z_1,\dots ,z_n),\sum _{k=0}^n{\mathbb{U}}_k(z_0,z_1,\dots ,z_k)\right\}\subseteq \{0,1\}.$$

$$({\mathbb{R}}^d{)}^{n+1}\ni (z_0,z_1,\dots ,z_n)\mapsto {\mathbb{U}}_n(z_0,z_1,\dots ,z_n)\in \{0,1\}$$

$$\forall k\in \{0,1,\dots ,n\}:{\mathbb{U}}_k(z_0,z_1,\dots ,z_k)={\mathbb{V}}_k(z_0,z_1,\dots ,z_k).$$

$$\sum_{k=0}^{N}\mathbb{V}_{k}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{k}(\omega)\big{)}=\sum_{k=0}^{N}\mathbbm{1}_{\{\tau=k\}}(\omega)=1.$$

$$\mathbb{U}_{k}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{k}(\omega)\big{)}=\mathbb{V}_{k}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{k}(\omega)\big{)}=\mathbbm{1}_{\{\tau=k\}}(\omega).$$

$$\tau(\omega)=\sum_{n=0}^{N}n\,\mathbb{U}_{n}\big{(}\mathcal{Z}_{0}(\omega),\mathcal{Z}_{1}(\omega),\dots,\mathcal{Z}_{n}(\omega)\big{)}.$$

$$\tau =\sum _{n=0}^{N}n{\mathbb{U}}_{n,\tau }({\mathcal{Z}}_0,{\mathcal{Z}}_1,\dots ,{\mathcal{Z}}_n)$$

$$U_{n,\theta }(z_0,z_1,\dots ,z_n)=\max \left\{u_{n,\theta }(z_n),n+1-N\right\}\left[1-\sum _{k=0}^{n-1}{U}_{k,\theta }(z_0,z_1,\dots ,z_k)\right]$$

$$\sum _{n=0}^N{U}_{n,\theta }(z_0,z_1,\dots ,z_n)=1.$$

Full text

[TABLE]

1. 1 Introduction

2. 2 Main ideas of the proposed algorithm

3. 2.1 The American option pricing problem

4. 2.2 Temporal discretisation

5. 2.3 Factorisation lemma for stopping times

6. 2.4 Neural network architectures for stopping times

7. 2.5 Formulation of the objective function

8. 2.6 Stochastic gradient ascent optimisation algorithms

9. 2.7 Price and optimal exercise time for American-style options

10. 3 Details of the proposed algorithm

11. 3.1 Formulation of the proposed algorithm in a special case

12. 3.2 Formulation of the proposed algorithm in the general case

13. 3.3 Comments on the proposed algorithm

14. 4 Numerical examples of pricing American-style derivatives

15. 4.1 Theoretical considerations

16. 4.1.1 Option prices in the Black–Scholes model

17. 4.1.2 Approximating American options with Bermudan options

18. 4.2 Setting

19. 4.3 Examples with known one-dimensional representation

20. 4.3.1 Optimal stopping of a Brownian motion

21. 4.3.1.1 A Bermudan three-exercise put-type example

22. 4.3.1.2 An American put-type example

23. 4.3.2 Geometric average-type options

24. 4.3.2.1 An American geometric average put-type example

25. 4.3.2.2 An American geometric average call-type example

26. 4.3.2.3 Another American geometric average call-type example

27. 4.4 Examples without known one-dimensional representation

28. 4.4.1 Max-call options

29. 4.4.1.1 A Bermudan max-call standard benchmark example

30. 4.4.1.2 A high-dimensional Bermudan max-call benchmark example

31. 4.4.1.3 Another Bermudan max-call example

32. 4.4.2 A strangle spread basket option

33. 4.4.3 A put basket option in Dupire’s local volatility model

34. 4.4.4 A path-dependent financial derivative

35. 5 Python source codes

36. 5.1 A Python source code for the algorithm

37. 5.2 A Python source code associated to Subsection 4.4.1.1

38. 5.3 A Python source code associated to Subsection 4.4.1.3

39. 5.4 A Python source code associated to Subsection 4.4.2

40. 5.5 A Python source code associated to Subsection 4.4.3

41. 5.6 A Python source code associated to Subsection 4.4.4

1 Introduction

Nowadays many financial derivatives, such as American or Bermudan options, are of early exercise type. Contrary to European options, the holder of such an option has the right to exercise before the time of maturity. In models from mathematical finance for the appropriate pricing of early exercise options, this aspect gives rise to optimal stopping problems. The dimension of such optimal stopping problems can often be quite high since it corresponds to the number of underlying assets. Due to the curse of dimensionality (cf. Bellman [11]), high-dimensional optimal stopping problems are, however, notoriously difficult to solve. Such optimal stopping problems can in nearly all cases not be solved explicitly and it is an active topic of research to design and analyse approximation methods which are capable of approximately solving possibly high-dimensional optimal stopping problems. Many different approaches for numerically solving optimal stopping problems and, in particular, American and Bermudan option pricing problems have been studied in the literature; cf., e.g., [89, 35, 5, 28, 1, 25, 90, 2, 27, 75, 91, 80, 83, 4, 43, 3, 23, 26, 44, 51, 38, 41, 46, 19, 68, 79, 29, 39, 55, 20, 22, 24, 63, 67, 71, 76, 15, 64, 66, 81, 12, 13, 70, 36, 54, 65, 14, 18, 82, 31, 49, 16, 17, 57, 73, 69, 21, 32, 85, 8, 9, 47, 86, 92, 7, 42, 48, 72, 10, 74, 30]. For example, such approaches include approximating the Snell envelope or continuation values (cf., e.g., [89, 5, 28, 75]), computing optimal exercise boundaries (cf., e.g., [2]), and dual methods (cf., e.g., [80, 51]). Whereas in [51, 66] artificial neural networks with one hidden layer were employed to approximate continuation values, more recently numerical approximation methods for American and Bermudan option pricing that are based on deep learning were introduced; cf., e.g., [86, 85, 9, 42, 10, 72, 30]. More precisely, in [86, 85] deep neural networks are used to approximately solve the corresponding obstacle partial differential equation problem, in [9] the corresponding optimal stopping problem is tackled directly with a deep learning-based algorithm, [42] applies an extension of the deep BSDE solver from [50, 37] to the corresponding reflected backward stochastic differential equation problem, [30] suggests a different deep learning-based algorithm that relies on discretising backward stochastic differential equations, and in [10, 72] deep neural network-based variants of the classical algorithm introduced by Longstaff & Schwartz [75] are examined.

In this work, we propose an algorithm for solving general possibly high-dimensional optimal stopping problems; cf. Framework 3.2 in Subsection 3.2. In spirit, it is similar to the algorithm introduced in [9]. The proposed algorithm is based on deep learning and computes both approximations of an optimal stopping strategy and the optimal expected pay-off associated with the considered optimal stopping problem. In the context of pricing early exercise options, these correspond to approximations of an optimal exercise strategy and the price of the considered option, respectively. The derivation and implementation of the proposed algorithm consist of essentially the following three steps.

1. (I)

A neural network architecture for, in an appropriate sense, ‘randomised’ stopping times (cf. (31) in Subsection 2.4) is established in such a way that varying the neural network parameters leads to different randomised stopping times being expressed. This neural network architecture is used to replace the supremum of the expected pay-off over suitable stopping times (which constitutes the generic optimal stopping problem) by the supremum of a suitable objective function over neural network parameters (cf. (2.5)–(39) in Subsection 2.5). 2. (II)

A stochastic gradient ascent-type optimisation algorithm is employed to compute neural network parameters that approximately maximise the objective function (cf. Subsection 2.6). 3. (III)

From these neural network parameters and the corresponding randomised stopping time, a true stopping time is constructed which serves as the approximation of an optimal stopping strategy (cf. (44) and (46) in Subsection 2.7). In addition, an approximation of the optimal expected pay-off is obtained by computing a Monte Carlo approximation of the expected pay-off under this approximately optimal stopping strategy (cf. (45) in Subsection 2.7).

It follows from (III) that the proposed algorithm computes a low-biased approximation of the optimal expected pay-off (cf. (48) in Subsection 2.7). Yet, a large number of numerical experiments where a reference value is available (cf. Section 4) show that the bias appears to become small quickly during training and that a very satisfying accuracy can be achieved in short computation time, even in high dimensions (cf. the end of this introduction below for a brief overview of the numerical computations that were performed). Moreover, in (I) we resort to randomised stopping times in order to circumvent the discrete nature of stopping times that attain only finitely many different values. As a result, it is possible in (II) to tackle the arising optimisation problem with a stochastic gradient ascent-type algorithm. Furthermore, while the focus in this article lies on American and Bermudan option pricing, the proposed algorithm can also be applied to optimal stopping problems that arise in other areas where the underlying stochastic process can be efficiently simulated. Apart from this, we only rely on the assumption that the stochastic process to be optimally stopped is a Markov process (cf. Subsection 2.4). But this assumption is no substantial restriction since, on the one hand, it is automatically fulfilled in many relevant problems and, on the other hand, a discrete stochastic process that is not a Markov process can be replaced by a Markov process of higher dimension that aggregates all necessary information (cf., e.g., [9, Subsection 4.3] and, e.g., Subsection 4.4.4).

Next we compare our algorithm to the one introduced in [9]. The latter splits the original problem into smaller optimal stopping problems at each time step where stopping is permitted and decides to stop at that point in time or later (cf. [9, (4) in Subsection 2.1]). Starting at maturity, these auxiliary problems are solved recursively backwards until the initial time is reached. Thereby, in every new step, neural network parameters are learned for an objective function that depends, in particular, on the parameters found in the previous steps (cf. [9, Subsection 2.3]). In contrast, in (I) a single objective function is designed. This objective function allows to search in (II) for neural network parameters that maximise the expected pay-off simultaneously over (randomised) stopping times which may decide to stop at any of the admissible points in time. Therefore, the algorithm proposed here does not rely on a recursion over the different time points. In addition, the construction of the final approximation of an optimal stopping strategy in (III) differs from a corresponding construction in [9]. We refer to Subsection 4.3.2.1 for a comparison between the two algorithms with respect to performance.

The remainder of this article is organised as follows. In Section 2, we present the main ideas from which the proposed algorithm is derived. More specifically, in Subsection 2.1 we illustrate how an optimal stopping problem in the context of American option pricing is typically formulated. Thereafter, a replacement of this continuous-time problem by a corresponding discrete time optimal stopping problem is discussed by means of an example in Subsection 2.2. Subsection 2.3 is devoted to the statement and proof of an elementary but crucial result about factorising general discrete stopping times in terms of compositions of measurable functions (cf. Lemma 2.2), which lies at the heart of the neural network architecture we propose in Subsection 2.4 to approximate general discrete stopping times. This construction, in turn, is exploited in Subsection 2.5 to transform the discrete time optimal stopping problem from Subsection 2.2 into the search of a maximum of a suitable objective function (cf. (I) above). In Subsection 2.6, we suggest to employ stochastic gradient ascent-type optimisation algorithms to find approximate maximum points of the objective function (cf. (II) above). As a last step, we explain in Subsection 2.7 how we calculate final approximations of both the American option price and an optimal exercise strategy (cf. (III) above). In Section 3, we introduce the proposed algorithm in a concise way, first for a special case for the sake of clarity (cf. Subsection 3.1) and second in more generality so that, in particular, a rigorous description of our implementations is fully covered (cf. Subsections 3.2–3.3). Following this, in Section 4 first a few theoretical results are presented (cf. Subsection 4.1), which are used to design numerical example problems and to provide reference values. Thereafter, we describe in detail a large number of example problems, on which our proposed algorithm was tested, and present numerical results for each of these problems. In particular, the examples include the optimal stopping of Brownian motions (cf. Subsection 4.3.1), the pricing of certain exotic American geometric average put and call-type options (cf. Subsection 4.3.2), the pricing of Bermudan max-call options in up to 5000 dimensions (cf. Subsection 4.4.1), the pricing of an American strangle spread basket option in five dimensions (cf. Subsection 4.4.2), the pricing of an American put basket option in Dupire’s local volatility model in five dimensions (cf. Subsection 4.4.3), and the pricing of an exotic path-dependent financial derivative of a single underlying, which is modelled as a 100-dimensional optimal stopping problem (cf. Subsection 4.4.4). The numerical results for the examples in Subsections 4.3.1.1, 4.3.1.2, 4.3.2.1, 4.3.2.2, 4.3.2.3, and 4.4.1.3 are compared to calculated reference values that can be easily obtained due to the specific design of the considered optimal stopping problem. Moreover, the examples in Subsections 4.3.2.2, 4.4.1.1, 4.4.1.3, 4.4.2, 4.4.3, and 4.4.4 are taken from the literature, and our corresponding numerical results are compared to reference values from the literature (where available). Finally, Python source codes employed to generate the numerical results in Subsections 4.4.1.1, 4.4.1.3, 4.4.2, 4.4.3, and 4.4.4 can be found in Section 5.

2 Main ideas of the proposed algorithm

In this section, we outline the main ideas that lead to the formulation of the proposed algorithm in Subsections 3.1–3.2 by considering the example of pricing an American option. The proposed algorithm in Framework 3.2 in Subsection 3.2 is, however, general enough to also be applied to optimal stopping problems where there are no specific assumptions on the dynamics of the underlying stochastic process, as long as it can be cheaply simulated (cf. Subsection 3.3). Furthermore, often in practice and, in particular, in the case of Bermudan option pricing (cf. many of the examples in Section 4), the optimal stopping problem of interest is not a continuous-time problem but is already formulated in discrete time. In such a situation, there is no need for a time discretisation, as described in Subsection 2.2 below, and the proposed algorithm in Framework 3.2 can be applied directly.

2.1 The American option pricing problem

Let $T\in(0,\infty)$, $d\in\mathbb{N}=\{1,2,3,\ldots\}$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space with a filtration $\mathscr{F}=(\mathscr{F}_{t})_{t\in[0,T]}$ that satisfies the usual conditions (cf., e.g., [59, Definition 2.25 in Section 1.2]), let $\xi\colon\Omega\to\mathbb{R}^{d}$ be an $\mathscr{F}_{0}$/$\mathcal{B}(\mathbb{R}^{d})$-measurable function which satisfies for all $p\in(0,\infty)$ that \mathbb{E}\bigl{[}\lVert\xi\rVert_{\mathbb{R}^{d}}^{p}\bigr{]}<\infty, let $W\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be a standard $(\Omega,\mathcal{F},\mathbb{P},\mathscr{F})$-Brownian motion with continuous sample paths, let $\mu\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ be Lipschitz continuous functions, let $X\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be an $\mathscr{F}$-adapted continuous solution process of the stochastic differential equation

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, and let $g\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ be a continuous and at most polynomially growing function. We think of $X$ as a model for the price processes of $d$ underlyings (say, $d$ stock prices) under the risk-neutral pricing measure $\mathbb{P}$ (cf., e.g., Kallsen [58]) and we are interested in approximatively pricing the American option on the process $(X_{t})_{t\in[0,T]}$ with the discounted pay-off function $g\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$, that is, we intend to compute the real number

[TABLE]

In addition to the price of the American option in the model (1), there is also a high demand from the financial engineering industry to compute an approximately optimal exercise strategy, that is, to compute a stopping time which approximately reaches the supremum in (2).

In a very simple example of (1)–(2), we can think of an American put option in the one-dimensional Black–Scholes model, in which there are an interest rate $r\in\mathbb{R}$, a dividend yield $\delta\in[0,\infty)$, a volatility $\beta\in(0,\infty)$, and a strike price $K\in(0,\infty)$ such that it holds for all $x\in\mathbb{R}$, $t\in[0,T]$ that $d=1$, $\mu(x)=(r-\delta)\,x$, $\sigma(x)=\beta\,x$, and $g(t,x)=e^{-rt}\max\{K-x,0\}$.

2.2 Temporal discretisation

To derive the proposed approximation algorithm, we first apply the Euler–Maruyama scheme to the stochastic differential equation (1) (cf. (5)–(6) below) and we employ a suitable time discretisation for the optimal stopping problem (2). For this let $N\in\mathbb{N}$ be a natural number and let $t_{0},t_{1},\dots,t_{N}\in[0,T]$ be real numbers with

[TABLE]

(such that the maximal mesh size $\max_{n\in\{0,1,\dots,N-1\}}(t_{n+1}-t_{n})$ is sufficiently small). Observe that (1) ensures that for all $n\in\{0,1,\dots,N-1\}$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Note that (4) suggests for every $n\in\{0,1,\dots,N-1\}$ that

[TABLE]

The approximation scheme associated with (5) is referred to as the Euler–Maruyama scheme in the literature (cf., e.g., Maruyama [77] and Kloeden & Platen [62]). More formally, let $\mathcal{X}=(\mathcal{X}^{(1)},\ldots,\mathcal{X}^{(d)})\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d}$ be the stochastic process which satisfies for all $n\in\{0,1,\dots,N-1\}$ that $\mathcal{X}_{0}=\xi$ and

[TABLE]

and let $\mathfrak{F}=(\mathfrak{F}_{n})_{n\in\{0,1,\dots,N\}}$ be the filtration generated by $\mathcal{X}$. Combining this with (5) suggests the approximation

[TABLE]

for the price (2) of the American option in Subsection 2.1. Below we employ, in particular, (7) to derive the proposed approximation algorithm.

For every $n\in\{0,1,\dots,N\}$ note that the current pay-off $g(t_{n},\mathcal{X}_{n})$ does not carry any information that is not already contained in $\mathcal{X}_{n}$. But typically, optimal exercise strategies are learned more efficiently if it is added as an additional feature. Therefore, we introduce the $(d+1)$-dimensional Markov process $\mathcal{Z}\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d+1}$ which satisfies for all $n\in\{0,1,\dots,N\}$ that \mathcal{Z}_{n}=\bigl{(}\mathcal{X}_{n}^{(1)},\dots,\mathcal{X}_{n}^{(d)},g(t_{n},\mathcal{X}_{n})\bigr{)}. Observe that $\mathcal{Z}$ and $\mathcal{X}$ generate the same filtration $\mathfrak{F}$.

2.3 Factorisation lemma for stopping times

The derivation of the proposed approximation algorithm is in parts based on an elementary reformulation of time-discrete stopping times (cf. the left-hand side of (7) above) in terms of measurable functions that appropriately characterise the behaviour of the stopping time; cf. (10) and (9) in Lemma 2.2 below. The proof of Lemma 2.2 employs the following well-known factorisation result, Lemma 2.1. Lemma 2.1 follows, e.g., from Klenke [61, Corollary 1.97].

Lemma 2.1** (Factorisation lemma).**

Let $(S,\mathcal{S})$ be a measurable space, let $\Omega$ be a set, let $B\in\mathcal{B}(\mathbb{R}\cup\{-\infty,\infty\})$, and let $X\colon\Omega\to S$ and $Y\colon\Omega\to B$ be functions. Then it holds that $Y$ is $\{X^{-1}(A)\colon A\in\mathcal{S}\}$/$\mathcal{B}(B)$-measurable if and only if there exists an $\mathcal{S}$/$\mathcal{B}(B)$-measurable function $f\colon S\to B$ such that

[TABLE]

We are now ready to present the above-mentioned Lemma 2.2. This elementary lemma is a consequence of Lemma 2.1 above.

Lemma 2.2** (Factorisation lemma for stopping times).**

Let $d,N\in\mathbb{N}$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, let $\mathcal{Z}\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d+1}$ be a stochastic process, and let $\mathfrak{F}=(\mathfrak{F}_{n})_{n\in\{0,1,\dots,N\}}$ be the filtration generated by $\mathcal{Z}$. Then

1. (i)

for all Borel measurable functions $\mathbb{U}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, with $\forall\,z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}\colon\sum_{n=0}^{N}\mathbb{U}_{n}(z_{0},z_{1},\dots,z_{n})=1$ it holds that the function

[TABLE]

is an $\mathfrak{F}$-stopping time and 2. (ii)

for every $\mathfrak{F}$-stopping time $\tau\colon\Omega\to\{0,1,\dots,N\}$ there exist Borel measurable functions $\mathbb{U}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, which satisfy $\forall\,z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}\colon\sum_{n=0}^{N}\mathbb{U}_{n}(z_{0},z_{1},\dots,z_{n})=1$ and

[TABLE]

Proof of Lemma 2.2.

Note that for all Borel measurable functions $\mathbb{U}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, with $\forall\,z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}\colon\sum_{n=0}^{N}\mathbb{U}_{n}(z_{0},z_{1},\dots,z_{n})=1$ and all $k\in\{0,1,\dots,N\}$ it holds that

[TABLE]

This establishes (i). It thus remains to prove (ii). For this let $\tau\colon\Omega\to\{0,1,\dots,N\}$ be an $\mathfrak{F}$-stopping time. Observe that for every function $\varrho\colon\Omega\to\{0,1,\dots,N\}$ and every $\omega\in\Omega$ it holds that

[TABLE]

Next note that for every $n\in\{0,1,\dots,N\}$ it holds that the function

[TABLE]

is $\mathfrak{F}_{n}$/$\mathcal{B}(\{0,1\})$-measurable. This and the fact that

[TABLE]

ensure that for every $n\in\{0,1,\dots,N\}$ it holds that the function

[TABLE]

is $\sigma_{\Omega}((\mathcal{Z}_{0},\mathcal{Z}_{1},\dots,\mathcal{Z}_{n}))$/$\mathcal{B}(\{0,1\})$-measurable. Lemma 2.1 hence demonstrates that there exist Borel measurable functions $\mathbb{V}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, which satisfy for all $n\in\{0,1,\dots,N\}$, $\omega\in\Omega$ that

[TABLE]

Next let $\mathbb{U}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\mathbb{R}$, $n\in\{0,1,\dots,N\}$, be the functions which satisfy for all $n\in\{0,1,\dots,N\}$, $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ that

[TABLE]

Observe that (17), in particular, ensures that for all $z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}$ it holds that

[TABLE]

Hence, we obtain that for all $z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}$ it holds that

[TABLE]

In addition, note that (17) assures that for all $z_{0}\in\mathbb{R}^{d+1}$ it holds that

[TABLE]

Induction, the fact that

[TABLE]

and (17) hence demonstrate that for all $n\in\{0,1,\dots,N\}$, $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ it holds that

[TABLE]

Moreover, note that (17), induction, and the fact that the functions $\mathbb{V}_{n}\colon(\mathbb{R}^{d})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, are Borel measurable ensure that for every $n\in\{0,1,\dots,N\}$ it holds that the function

[TABLE]

is also Borel measurable. In the next step, we observe that (20), (17), (21), and induction assure that for all $n\in\{0,1,\dots,N\}$, $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ with $n+1-N\leq\sum_{k=0}^{n}\mathbb{V}_{k}(z_{0},z_{1},\dots,z_{k})\leq 1$ it holds that

[TABLE]

In addition, note that (16) shows that for all $\omega\in\Omega$ it holds that

[TABLE]

This, (24), and again (16) imply that for all $k\in\{0,1,\dots,N\}$, $\omega\in\Omega$ it holds that

[TABLE]

Equation (12) hence proves that for all $\omega\in\Omega$ it holds that

[TABLE]

Combining this with (19) and (23) establishes (ii). The proof of Lemma 2.2 is thus complete. ∎

2.4 Neural network architectures for stopping times

In the next step, we employ multilayer neural network approximations of the functions $\mathbb{U}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, in the factorisation lemma, Lemma 2.2 above. In the following, we refer to these functions as ‘stopping time factors’. Consider again the setting in Subsections 2.1–2.2, for every $\mathfrak{F}$-stopping time $\tau\colon\Omega\to\{0,1,\dots,N\}$ let $\mathbb{U}_{n,\tau}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$, $n\in\{0,1,\dots,N\}$, be Borel measurable functions which satisfy $\forall\,z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}\colon\sum_{n=0}^{N}\mathbb{U}_{n,\tau}(z_{0},z_{1},\dots,z_{n})=1$ and

[TABLE]

(cf. (ii) of Lemma 2.2), let $\nu\in\mathbb{N}$ be a sufficiently large natural number, and for every $n\in\{0,1,\dots,N\}$, $\theta\in\mathbb{R}^{\nu}$ let $u_{n,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$ and $U_{n,\theta}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$ be Borel measurable functions which satisfy for all $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ that

[TABLE]

(cf. (17) above). We think of $\nu\in\mathbb{N}$ as the number of parameters in the employed artificial neural networks (for instance, in the case of (32)–(33) below we have $\nu=N(l^{2}+l(d+4)+1)$, where $l$ is the number of neurons in each of the two hidden layers), for every $n\in\{0,1,\dots,N\}$, $\theta\in\mathbb{R}^{\nu}$ we think of $u_{n,\theta}\colon\mathbb{R}^{d+1}\rightarrow(0,1)$ as an appropriate artificial neural network (cf. the next paragraph for further details), and for every appropriate $\mathfrak{F}$-stopping time $\tau\colon\Omega\to\{0,1,\dots,N\}$ and every $n\in\{0,1,\dots,N\}$ we think of the function $U_{n,\theta}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$ for suitable $\theta\in\mathbb{R}^{\nu}$ as an appropriate approximation of the stopping time factor $\mathbb{U}_{n,\tau}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$. In addition, it shall be noted that for all $n\in\{0,1,\dots,N\}$, in comparison with the function $\mathbb{V}_{n}\colon(\mathbb{R}^{d+1})^{n+1}\to\{0,1\}$ in (17) above, the functions $u_{n,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, in (29) are defined only on $\mathbb{R}^{d+1}$ instead of $(\mathbb{R}^{d+1})^{n+1}$ and, therefore, only depend on $z_{n}\in\mathbb{R}^{d+1}$ instead of the whole vector $(z_{0},z_{1},\ldots,z_{n})\in(\mathbb{R}^{d+1})^{n+1}$. Even though this constitutes a significant simplification, [9, Theorem 1 and Remark 2 in Subsection 2.1] suggest that, due to the fact that $\mathcal{Z}\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d+1}$ is a Markov process, the approximate stopping time factors $U_{n,\theta}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, $n\in\{0,1,\dots,N\}$, still possess enough flexibility to represent an optimal stopping time for the discrete stopping problem corresponding to the left side of (7). Furthermore, observe that for all $\theta\in\mathbb{R}^{\nu}$, $z_{0},z_{1},\dots,z_{N}\in\mathbb{R}^{d+1}$ it holds that

[TABLE]

Because of this, for every $\theta\in\mathbb{R}^{\nu}$ the stochastic process

[TABLE]

can also be viewed as an appropriate sense ‘randomised stopping time’ (cf., e.g., [88, Definition 1 in Subsection 3.1] and, e.g., [40, Section 1.1]).

We suggest to choose the functions $u_{n,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, $n\in\{0,1,\dots,N-1\}$, as multilayer feedforward neural networks (cf. [9, Corollary 5 in Subsection 2.2] and, e.g., [33, 52, 6]). For example, for every $k\in\mathbb{N}$ let $\mathcal{L}_{k}\colon\mathbb{R}^{k}\to\mathbb{R}^{k}$ be the function which satisfies for all $x=(x_{1},\dots,x_{k})\in\mathbb{R}^{k}$ that

[TABLE]

for every $\theta=(\theta_{1},\dots,\theta_{\nu})\in\mathbb{R}^{\nu}$, $v\in\mathbb{N}_{0}=\{0,1,2,\ldots\}$, $k,j\in\mathbb{N}$ with $v+k(j+1)\leq\nu$ let $A^{\theta,v}_{k,j}\colon\mathbb{R}^{j}\to\mathbb{R}^{k}$ be the affine linear function which satisfies for all $x=(x_{1},\dots,x_{j})\in\mathbb{R}^{j}$ that

[TABLE]

let $l\in\mathbb{N}$, and assume for all $n\in\{0,1,\dots,N-1\}$, $\theta\in\mathbb{R}^{\nu}$ that $\nu\geq N(l^{2}+l(d+4)+1)$ and

[TABLE]

The functions in (34) provide artificial neural networks with $4$ layers ($1$ input layer with $d+1$ neurons, $2$ hidden layers with $l$ neurons each, and $1$ output layer with $1$ neuron) and the multidimensional version of the standard logistic function $\mathbb{R}\ni x\mapsto\nicefrac{{\exp(x)}}{{(\exp(x)+1)}}\in(0,1)$ (cf. (32) above) as activation functions. In our numerical simulations in Section 4, we use this type of activation function only just in front of the output layer and we employ instead the multidimensional version of the rectifier function $\mathbb{R}\ni x\mapsto\max\{x,0\}\in[0,\infty)$ as activation functions just in front of the hidden layers. But in order to keep the illustration here as short as possible, we only employ the multidimensional version of the standard logistic function as activation functions in (32)–(34) above. Furthermore, note that in contrast to the choice of the functions $u_{n,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, $n\in\{0,1,\dots,N-1\}$, the choice of the functions $u_{N,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, has no influence on the approximate stopping time factors $U_{n,\theta}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$, $\theta\in\mathbb{R}^{\nu}$, $n\in\{0,1,\dots,N\}$ (cf. (29) above).

2.5 Formulation of the objective function

Recall that we intend to compute the real number

[TABLE]

as an approximation of the American option price (2) (cf. (7) in Subsection 2.2). By employing neural network architectures for stopping times (cf. Subsection 2.4 above), we next propose to replace the search over all $\mathfrak{F}$-stopping times for finding the supremum in (35) by a search over the artificial neural network parameters $\theta\in\mathbb{R}^{\nu}$ (cf. (2.5) below). For this, observe that (28) implies for all $\mathfrak{F}$-stopping times $\tau\colon\Omega\to\{0,1,\dots,N\}$ and all $n\in\{0,1,\dots,N\}$ that

[TABLE]

Therefore, for all $\mathfrak{F}$-stopping times $\tau\colon\Omega\to\{0,1,\dots,N\}$ it holds that

[TABLE]

Combining this with (i) of Lemma 2.2 and (30) inspires the approximation

[TABLE]

In view of this, our numerical solution for approximatively computing (35) consists of trying to find an approximate maximiser of the objective function

[TABLE]

2.6 Stochastic gradient ascent optimisation algorithms

Local/global maxima of the objective function (39) can be approximately reached by maximising the expectation of the random objective function

[TABLE]

by means of a stochastic gradient ascent-type optimisation algorithm. This yields a sequence of random parameter vectors along which we expect the objective function (39) to increase. More formally, applying under suitable hypotheses stochastic gradient ascent-type optimisation algorithms to (39) results in random approximations

[TABLE]

for $m\in\{0,1,2,\dots\}$ of the local/global maximum points of the objective function (39), where $m\in\{0,1,2,\dots\}$ is the number of steps of the employed stochastic gradient ascent-type optimisation algorithm.

2.7 Price and optimal exercise time for American-style options

The approximation algorithm sketched in Subsection 2.6 above allows us to approximatively compute both the price and an optimal exercise strategy for the American option (cf. Subsection 2.1). Let $M\in\mathbb{N}$ and consider a realisation $\widehat{\Theta}_{M}\in\mathbb{R}^{\nu}$ of the random variable $\Theta_{M}\colon\Omega\to\mathbb{R}^{\nu}$. Then for sufficiently large $N,\nu,M\in\mathbb{N}$ a candidate for a suitable approximation of the American option price is the real number

[TABLE]

and a candidate for a suitable approximation of an optimal exercise strategy for the American option is the function

[TABLE]

Note, however, that in general the function (43) does not take values in $\{0,1,\ldots,N\}$ and hence is not a proper stopping time. Similarly, note that in general it is not clear whether there exists an exercise strategy such that the number (42) is equal to the expected discounted pay-off under this exercise strategy. For these reasons, we suggest other candidates for suitable approximations of the price and an optimal exercise strategy for the American option. More specifically, for every $\theta\in\mathbb{R}^{\nu}$ let $\tau_{\theta}\colon\Omega\to\{0,1,\dots,N\}$ be the $\mathfrak{F}$-stopping time given by

[TABLE]

(cf. (30) above). Then for sufficiently large $N,\nu,M\in\mathbb{N}$ we use a suitable Monte Carlo approximation of the real number

[TABLE]

as a suitable implementable approximation of the price of the American option (cf. (2) in Subsection 2.1 above and (59) in Subsection 3.1 below) and we use the random variable

[TABLE]

as a suitable implementable approximation of an optimal exercise strategy for the American option. Note that (30) ensures that

[TABLE]

This shows that the exercise strategy $\tau_{\widehat{\Theta}_{M}}\colon\Omega\to\{0,1,\dots,N\}$ exercises at the first time index $n\in\{0,1,\dots,N\}$ for which the approximate stopping time factor associated with the mesh point $t_{n}$ is at least as large as the combined approximate stopping time factors associated with all later mesh points. Finally, observe that

[TABLE]

This implies that Monte Carlo approximations of the number (45) typically are low-biased approximations of the American option price (2).

3 Details of the proposed algorithm

3.1 Formulation of the proposed algorithm in a special case

In this subsection, we describe the proposed algorithm in the specific situation where the objective is to solve the American option pricing problem described in Subsection 2.1, where batch normalisation (cf. Ioffe & Szegedy [53]) is not employed in the proposed algorithm, and where the plain vanilla stochastic gradient ascent approximation method with a constant learning rate $\gamma\in(0,\infty)$ and without mini-batches is the employed stochastic approximation algorithm. The general framework, which includes the setting in this subsection as a special case, can be found in Subsection 3.2 below.

Framework 3.1** (Specific case).**

Let $T,\gamma\in(0,\infty)$, $d,N,l\in\mathbb{N}$, $\nu=N(l^{2}+l(d+4)+1)$, let $\mu\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$, $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$, and $g\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ be Borel measurable functions, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, let $\xi^{m}\colon\Omega\to\mathbb{R}^{d}$, $m\in\mathbb{N}$, be independent random variables, let $W^{m}\colon[0,T]\times\Omega\to\mathbb{R}^{d}$, $m\in\mathbb{N}$, be independent $\mathbb{P}$-standard Brownian motions with continuous sample paths, assume that $(\xi^{m})_{m\in\mathbb{N}}$ and $(W^{m})_{m\in\mathbb{N}}$ are independent, let $t_{0},t_{1},\dots,t_{N}\in[0,T]$ be real numbers with $0=t_{0}<t_{1}<\ldots<t_{N}=T$, let $\mathcal{X}^{m}=(\mathcal{X}^{m,(1)},\ldots,\mathcal{X}^{m,(d)})\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d}$, $m\in\mathbb{N}$, and $\mathcal{Z}^{m}\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d+1},$ $m\in\mathbb{N}$, be the stochastic processes which satisfy for all $m\in\mathbb{N}$, $n\in\{0,1,\dots,N-1\}$, $\mathfrak{n}\in\{0,1,\dots,N\}$ that $\mathcal{X}^{m}_{0}=\xi^{m}$,

[TABLE]

and

[TABLE]

for every $k\in\mathbb{N}$ let $\mathcal{L}_{k}\colon\mathbb{R}^{k}\to\mathbb{R}^{k}$ be the function which satisfies for all $x=(x_{1},\dots,x_{k})\in\mathbb{R}^{k}$ that

[TABLE]

for every $\theta=(\theta_{1},\dots,\theta_{\nu})\in\mathbb{R}^{\nu}$, $v\in\mathbb{N}_{0}$, $k,j\in\mathbb{N}$ with $v+k(j+1)\leq\nu$ let $A^{\theta,v}_{k,j}\colon\mathbb{R}^{j}\to\mathbb{R}^{k}$ be the function which satisfies for all $x=(x_{1},\dots,x_{j})\in\mathbb{R}^{j}$ that

[TABLE]

for every $\theta\in\mathbb{R}^{\nu}$ let $u_{n,\theta}\colon\mathbb{R}^{d+1}\to(0,1)$, $n\in\{0,1,\dots,N\}$, be functions which satisfy for all $n\in\{0,1,\dots,N-1\}$ that

[TABLE]

for every $n\in\{0,1,\dots,N\}$, $\theta\in\mathbb{R}^{\nu}$ let $U_{n,\theta}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$ be the function which satisfies for all $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ that

[TABLE]

for every $m\in\mathbb{N}$ let $\phi^{m}\colon\mathbb{R}^{\nu}\times\Omega\to\mathbb{R}$ be the function which satisfies for all $\theta\in\mathbb{R}^{\nu}$, $\omega\in\Omega$ that

[TABLE]

for every $m\in\mathbb{N}$ let $\Phi^{m}\colon\mathbb{R}^{\nu}\times\Omega\to\mathbb{R}^{\nu}$ be the function which satisfies for all $\theta\in\mathbb{R}^{\nu}$, $\omega\in\Omega$ that

[TABLE]

let $\Theta\colon\mathbb{N}_{0}\times\Omega\to\mathbb{R}^{\nu}$ be a stochastic process which satisfies for all $m\in\mathbb{N}$ that

[TABLE]

and for every $j\in\mathbb{N}$, $\theta\in\mathbb{R}^{\nu}$ let $\tau_{j,\theta}\colon\Omega\to\{0,1,\dots,N\}$ be the random variable given by

[TABLE]

Consider the setting of Framework 3.1, assume that $\mu$ and $\sigma$ are globally Lipschitz continuous, and assume that $g$ is continuous and at most polynomially growing. In the case of sufficiently large $N,M,J\in\mathbb{N}$ and sufficiently small $\gamma\in(0,\infty)$, we then think of the random number

[TABLE]

as an approximation of the price of the American option with the discounted pay-off function $g$ and for every $j\in\mathbb{N}$ we think of the random variable

[TABLE]

as an approximation of an optimal exercise strategy associated with the underlying time-discrete path $(\mathcal{X}^{M+j}_{n})_{n\in\{0,1,\dots,N\}}$ (cf. Subsection 2.1 above and Section 4 below).

3.2 Formulation of the proposed algorithm in the general case

In this subsection, we extend the framework in Subsection 3.1 above and describe the proposed algorithm in the general case.

Framework 3.2**.**

Let $T\in(0,\infty)$, $d,N,M,\nu,\varsigma,\varrho\in\mathbb{N}$, let $g\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ be a Borel measurable function, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, let $t_{0},t_{1},\dots,t_{N}\in[0,T]$ be real numbers with $0=t_{0}<t_{1}<\ldots<t_{N}=T$, let $\mathcal{X}^{m,j}=(\mathcal{X}^{m,j,(1)},\ldots,\mathcal{X}^{m,j,(d)})\colon\allowbreak\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d}$, $m\in\mathbb{N}_{0}$, $j\in\mathbb{N}$, be i.i.d. stochastic processes, for every $m\in\mathbb{N}_{0}$, $j\in\mathbb{N}$ let $\mathcal{Z}^{m,j}\colon\{0,1,\dots,N\}\times\Omega\to\mathbb{R}^{d+1}$ be the stochastic process which satisfies for all $n\in\{0,1,\dots,N\}$ that \mathcal{Z}^{m,j}_{n}=\big{(}\mathcal{X}^{m,j,(1)}_{n},\dots,\mathcal{X}^{m,j,(d)}_{n},g(t_{n},\mathcal{X}^{m,j})\big{)}, for every $n\in\{0,1,\dots,N\}$, $\theta\in\mathbb{R}^{\nu}$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$ let $u_{n}^{\theta,\mathbf{s}}\colon\mathbb{R}^{d+1}\to(0,1)$ be a function, for every $n\in\{0,1,\dots,N\}$, $\theta\in\mathbb{R}^{\nu}$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$ let $U_{n}^{\theta,\mathbf{s}}\colon(\mathbb{R}^{d+1})^{n+1}\to(0,1)$ be the function which satisfies for all $z_{0},z_{1},\dots,z_{n}\in\mathbb{R}^{d+1}$ that

[TABLE]

let $(J_{m})_{m\in\mathbb{N}_{0}}\subseteq\mathbb{N}$ be a sequence, for every $m\in\mathbb{N}$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$ let $\phi^{m,\mathbf{s}}\colon\mathbb{R}^{\nu}\times\Omega\to\mathbb{R}$ be the function which satisfies for all $\theta\in\mathbb{R}^{\nu}$, $\omega\in\Omega$ that

[TABLE]

for every $m\in\mathbb{N}$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$ let $\Phi^{m,\mathbf{s}}\colon\mathbb{R}^{\nu}\times\Omega\to\mathbb{R}^{\nu}$ be a function which satisfies for all $\omega\in\Omega$, $\theta\in\{\eta\in\mathbb{R}^{\nu}\colon\phi^{m,\mathbf{s}}(\cdot,\omega)\colon\mathbb{R}^{\nu}\to\mathbb{R}\text{ is differentiable at }\eta\}$ that

[TABLE]

let $\mathcal{S}\colon\mathbb{R}^{\varsigma}\times\mathbb{R}^{\nu}\times(\mathbb{R}^{d+1})^{\{0,1,\ldots,N-1\}\times\mathbb{N}}\to\mathbb{R}^{\varsigma}$ be a function, for every $m\in\mathbb{N}$ let $\Psi_{m}\colon\mathbb{R}^{\varrho}\times\mathbb{R}^{\nu}\to\mathbb{R}^{\varrho}$ and $\psi_{m}\colon\mathbb{R}^{\varrho}\to\mathbb{R}^{\nu}$ be functions, let $\mathbb{S}\colon\mathbb{N}_{0}\times\Omega\to\mathbb{R}^{\varsigma}$, $\Xi\colon\mathbb{N}_{0}\times\Omega\to\mathbb{R}^{\varrho}$, and $\Theta\colon\mathbb{N}_{0}\times\Omega\to\mathbb{R}^{\nu}$ be stochastic processes which satisfy for all $m\in\mathbb{N}$ that

[TABLE]

[TABLE]

for every $j\in\mathbb{N}$, $\theta\in\mathbb{R}^{\nu}$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$ let $\tau^{j,\theta,\mathbf{s}}\colon\Omega\to\{0,1,\dots,N\}$ be the random variable given by

[TABLE]

and let $\mathcal{P}\colon\Omega\to\mathbb{R}$ be the random variable which satisfies for all $\omega\in\Omega$ that

[TABLE]

Consider the setting of Framework 3.2. Under suitable further assumptions, in the case of sufficiently large $N,M,\nu,J_{0}\in\mathbb{N}$, we think of the random number

[TABLE]

as an approximation of the price of the American option with the discounted pay-off function $g$ and for every $j\in\mathbb{N}$ we think of the random variable

[TABLE]

as an approximation of an optimal exercise strategy associated with the underlying time-discrete path $(\mathcal{X}^{0,j}_{n})_{n\in\{0,1,\dots,N\}}$ (cf. Subsection 2.1 above and Section 4 below).

3.3 Comments on the proposed algorithm

Note that the lack in Framework 3.2 of any assumptions on the dynamics of the stochastic process $(\mathcal{X}_{n}^{0,1})_{n\in\{0,1,\dots,N\}}$ allows us to approximatively compute the optimal pay-off as well as an optimal exercise strategy for very general optimal stopping problems where, in particular, the stochastic process under consideration is not necessarily related to the solution of a stochastic differential equation. We only require that $(\mathcal{X}_{n}^{0,1})_{n\in\{0,1,\dots,N\}}$ can be simulated efficiently and formally we still rely on the Markov assumption (cf. Subsection 2.4 above). In addition, observe that the choice of the functions $u_{N}^{\theta,\mathbf{s}}\colon\mathbb{R}^{d+1}\to(0,1)$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$, $\theta\in\mathbb{R}^{\nu}$, has no influence on the proposed algorithm (cf. (61)). Furthermore, the dynamics in (65) associated with the stochastic processes $(\Xi_{m})_{m\in\mathbb{N}_{0}}$ and $(\Theta_{m})_{m\in\mathbb{N}_{0}}$ allow us to incorporate different stochastic approximation algorithms such as

• •

plain vanilla stochastic gradient ascent with or without mini-batches (cf. (57) above) as well as

• •

adaptive moment estimation (Adam) with mini-batches (cf. Kingma & Ba [60] and (97)–(98) in Subsection 4.2 below)

into the algorithm in Subsection 3.2 (cf. E, Han, & Jentzen [37, Subsection 3.3]). The dynamics in (64) associated with the stochastic process $(\mathbb{S}_{m})_{m\in\mathbb{N}_{0}}$, in turn, allow us to incorporate batch normalisation (cf. Ioffe & Szegedy [53] and the beginning of Section 4 below) into the algorithm in Subsection 3.2. In that case, we think of $(\mathbb{S}_{m})_{m\in\mathbb{N}_{0}}$ as a bookkeeping process keeping track of approximatively calculated means and standard deviations as well as of the number of steps $m\in\mathbb{N}_{0}$ of the employed stochastic approximation algorithm.

4 Numerical examples of pricing American-style derivatives

In this section, we test the algorithm of Framework 3.2 on several examples of pricing American-style financial derivatives. In each of the examples below, we employ the general approximation algorithm of Framework 3.2 above in conjunction with the Adam optimiser (cf. Kingma & Ba [60]) with varying learning rates and with mini-batches (cf. Subsection 4.2 below for a precise description).

In the example in Subsection 4.4.4 below, the initial value $\mathcal{X}_{0}^{0,1}$ is random. Therefore, we use $N$ fully connected feedforward neural networks to model the functions $u_{0}^{\theta,\mathbf{s}},\ldots,u_{N-1}^{\theta,\mathbf{s}}\colon\allowbreak\mathbb{R}^{d+1}\to(0,1)$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$, $\theta\in\mathbb{R}^{\nu}$. However, in all the other examples, $\mathcal{X}_{0}^{0,1}$ is deterministic. So it is enough to learn $N-1$ networks describing the functions $u_{1}^{\theta,\mathbf{s}},\ldots,u_{N-1}^{\theta,\mathbf{s}}\colon\mathbb{R}^{d+1}\to(0,1)$, $\mathbf{s}\in\mathbb{R}^{\varsigma}$, $\theta\in\mathbb{R}^{\nu}$. Then it can be decided whether it is better to stop at time [math] or not by comparing the deterministic pay-off $g(0,\mathcal{X}_{0}^{0,1})$ to a standard Monte Carlo estimate of the expected pay-off generated by the stopping strategy given by $u_{0}^{\theta,\mathbf{s}}=0$ and the functions $u_{1}^{\theta,\mathbf{s}},\ldots,u_{N-1}^{\theta,\mathbf{s}}\colon\mathbb{R}^{d+1}\to(0,1)$; cf. [9, Remark 6 in Subsection 2.3].

The standard network architecture we use in this paper consists of a $(d+1)$-dimensional input layer, two $(d+50)$-dimensional hidden layers, and a one-dimensional output layer. As non-linear activation functions just in front of the hidden layers, we employ the multidimensional version of the rectifier function $\mathbb{R}\ni x\mapsto\max\{x,0\}\in[0,\infty)$, whereas just in front of the output layer we employ the standard logistic function $\mathbb{R}\ni x\mapsto\nicefrac{{\exp(x)}}{{(\exp(x)+1)}}\in(0,1)$. In addition, batch normalisation (cf. Ioffe & Szegedy [53]) is applied just before the first linear transformation, just before each of the non-linear activation functions in front of the hidden layers as well as just before the non-linear activation function in front of the output layer. We use Xavier initialisation (cf. Glorot & Bengio [45]) to initialise all weights in the neural networks.

Two hidden layers work well in all our examples. However, the examples in Subsection 4.3 have an underlying one-dimensional structure, and as a consequence, fewer hidden layers yield equally good results; see Tables 2–3 below. On the other hand, the examples in Subsection 4.4 are more complex. In particular, it can be seen from Table 11 that for the max-call option in Subsection 4.4.1.1, two hidden layers give better results than zero or one hidden layer, but more than two hidden layers do not improve the results.

All examples presented below were implemented in Python. The correspondingPython source codes (cf. Section 5) were run, unless stated otherwise (cf. Subsection 4.4.1.2 as well as the last sentence in Subsection 4.4.1.3 below), in single precision (float32) on a NVIDIA GeForce RTX 2080 Ti GPU. The underlying system consisted of an AMD Ryzen 9 3950X CPU with 64 GB DDR4 memory running Tensorflow 2.1 on Ubuntu 19.10. We would like to point out that no special emphasis was put on optimising computation speed. In many cases, some of the algorithm parameters could be adjusted in order to obtain similarly accurate results in shorter runtime.

4.1 Theoretical considerations

Before we present the optimal stopping problem examples on which we tested the algorithm of Framework 3.2 (cf. Subsections 4.3–4.4 below), we recall a few theoretical results, which are used to design some of these examples, determine reference values, and provide further insights.

4.1.1 Option prices in the Black–Scholes model

The elementary and well-known result in Lemma 4.1 below specifies the distributions of linear combinations of independent and identically distributed centred Gaussian random variables which take values in a separable normed $\mathbb{R}$-vector space.

Lemma 4.1**.**

Let $n\in\mathbb{N}$, $\gamma=(\gamma_{1},\ldots,\gamma_{n})\in\mathbb{R}^{n}$, let $(V,\left\|\cdot\right\|_{V})$ be a separable normed $\mathbb{R}$-vector space, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, and let $X_{i}\colon\Omega\to V$, $i\in\{1,\ldots,n\}$, be i.i.d. centred Gaussian random variables. Then it holds that

[TABLE]

Proof of Lemma 4.1.

Throughout this proof let $Y_{1},Y_{2}\colon\Omega\to V$ be the random variables given by $Y_{1}=\sum_{i=1}^{n}\gamma_{i}\,X_{i}$ and $Y_{2}=\|\gamma\|_{\mathbb{R}^{n}}\,X_{1}$. Note that for every continuous linear functional $\varphi\colon V\to\mathbb{R}$ it holds that $\varphi\circ X_{i}\colon\Omega\to\mathbb{R}$, $i\in\{1,\ldots,n\}$, are independent and identically distributed centred Gaussian random variables. This implies for all continuous linear functionals $\varphi\colon V\to\mathbb{R}$ that

[TABLE]

This and, e.g., Jentzen, Salimova, & Welti [56, Lemma 4.10] establish that $Y_{1}(\mathbb{P})_{\mathcal{B}(V)}=Y_{2}(\mathbb{P})_{\mathcal{B}(V)}$. The proof of Lemma 4.1 is thus complete. ∎

The next elementary and well-known corollary follows directly from Lemma 4.1.

Corollary 4.2**.**

Let $d\in\mathbb{N}$, $\gamma=(\gamma_{1},\ldots,\gamma_{d})\in\mathbb{R}^{d}$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, and let $W=(W^{(1)},\ldots,W^{(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be a $\mathbb{P}$-standard Brownian motion with continuous sample paths. Then it holds that

[TABLE]

The next elementary result, Proposition 4.3, states that the distribution of a product of multiple correlated geometric Brownian motions is equal to the distribution of a single particular geometric Brownian motion.

Proposition 4.3**.**

Let $T,\epsilon\in(0,\infty)$, $d\in\mathbb{N}$, $\mathfrak{S}=(\varsigma_{1},\ldots,\varsigma_{d})\in\mathbb{R}^{d\times d}$, $\xi=(\xi_{1},\ldots,\xi_{d}),\alpha=(\alpha_{1},\ldots,\alpha_{d}),\beta=(\beta_{1},\ldots,\beta_{d})\in\mathbb{R}^{d}$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, let $\mathscr{F}^{(i)}=\linebreak(\mathscr{F}_{t}^{(i)})_{t\in[0,T]}$, $i\in\{1,2\}$, be filtrations on $(\Omega,\mathcal{F},\mathbb{P})$ that satisfy the usual conditions, let $W=(W^{(1)},\ldots,W^{(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be a standard $(\Omega,\mathcal{F},\mathbb{P},\mathscr{F}^{(1)})$-Brownian motion with continuous sample paths, let $\textsc{w}\colon[0,T]\times\Omega\to\mathbb{R}$ be a standard $(\Omega,\mathcal{F},\mathbb{P},\mathscr{F}^{(2)})$-Brownian motion with continuous sample paths, let $\mu\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$, $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$, $\mathbf{P}\colon C([0,T],\mathbb{R}^{d})\to C([0,T],\mathbb{R})$, and $\mathbf{G}\colon C([0,T],\mathbb{R})\to C([0,T],\mathbb{R})$ be the functions which satisfy for all $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$, $u^{(1)}=(u_{s}^{(1)})_{s\in[0,T]},\ldots,u^{(d)}=(u_{s}^{(d)})_{s\in[0,T]}\in C([0,T],\mathbb{R})$, $t\in[0,T]$ that $\mu(x)=(\alpha_{1}x_{1},\ldots,\alpha_{d}x_{d})$, $\sigma(x)=\operatorname{diag}(\beta_{1}x_{1},\ldots,\beta_{d}x_{d})\,\mathfrak{S}^{*}$, (\mathbf{G}[u^{(1)}])_{t}=\exp\bigl{(}\epsilon\,\bigl{[}\sum_{i=1}^{d}\alpha_{i}-\nicefrac{{\|\beta_{i}\varsigma_{i}\|_{\mathbb{R}^{d}}^{2}}}{{2}}\bigr{]}t\mathop{+}\epsilon\,\|\mathfrak{S}\,\beta\|_{\mathbb{R}^{d}}\,u_{t}^{(1)}\bigr{)}\prod_{i=1}^{d}|\xi_{i}|^{\epsilon}, and (\mathbf{P}[(u^{(1)},\ldots,u^{(d)})])_{t}=\prod_{i=1}^{d}\bigl{|}u_{t}^{(i)}\bigr{|}^{\epsilon}, let $X=\linebreak(X^{(1)},\ldots,X^{(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be an $\mathscr{F}^{(1)}$-adapted stochastic process with continuous sample paths, let $Y\colon[0,T]\times\Omega\to\mathbb{R}$ be an $\mathscr{F}^{(2)}$-adapted stochastic process with continuous sample paths, and assume that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Then

1. (i)

for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE] 2. (ii)

it holds that $\mathbf{P}$ and $\mathbf{G}$ are continuous functions, and 3. (iii)

it holds that

[TABLE]

Proof of Proposition 4.3.

Throughout this proof let $\gamma=(\gamma_{1},\ldots,\gamma_{d})\in\mathbb{R}^{d}$ be the vector given by $\gamma=\mathfrak{S}\,\beta$, let $Z^{(i)}\colon[0,T]\times\Omega\to\mathbb{R}$, $i\in\{1,\ldots,d\}$, be the stochastic processes which satisfy for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ that $Z_{t}^{(i)}=[\alpha_{i}-\nicefrac{{\|\beta_{i}\varsigma_{i}\|_{\mathbb{R}^{d}}^{2}}}{{2}}]t+\beta_{i}\langle\varsigma_{i},W_{t}\rangle_{\mathbb{R}^{d}}$, and let $\tilde{\mathbf{G}}\colon C([0,T],\mathbb{R})\to C([0,T],\mathbb{R})$ be the function which satisfies for all $u=(u_{s})_{s\in[0,T]}\in C([0,T],\mathbb{R})$, $t\in[0,T]$ that

[TABLE]

Observe that for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

In addition, note that for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Itô’s formula hence shows that for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Combining this and (78) with, e.g., Da Prato & Zabczyk [34, (i) in Theorem 7.4] proves that for all $i\in\{1,\ldots,d\}$, $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

This establishes (i). In the next step note that (ii) is clear. It thus remains to prove (iii). For this observe that (i) establishes that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Continuity hence implies that it holds $\mathbb{P}$-a.s. that

[TABLE]

Moreover, note that (i) shows that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

This and continuity establish that it holds $\mathbb{P}$-a.s. that

[TABLE]

Furthermore, observe that Corollary 4.2 ensures that

[TABLE]

The fact that $\tilde{\mathbf{G}}\colon C([0,T],\mathbb{R})\to C([0,T],\mathbb{R})$ is a Borel measurable function, (ii), the fact that $\forall\,u\in C([0,T],\mathbb{R})\colon\tilde{\mathbf{G}}[\|\gamma\|_{\mathbb{R}^{d}}\,u]=\mathbf{G}[u]$, (83), and (85) hence demonstrate that

[TABLE]

The proof of Proposition 4.3 is thus complete. ∎

In the next result, Lemma 4.4, we recall the well-known formula for the price of a European call option on a single stock in the Black–Scholes model (cf., e.g., Øksendal [78, Corollary 12.3.8]).

Lemma 4.4**.**

Let $T,\xi,\sigma\in(0,\infty)$, $r,c\in\mathbb{R}$, let $\Phi\colon\mathbb{R}\to\mathbb{R}$ be the function which satisfies for all $x\in\mathbb{R}$ that $\Phi(x)=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}\,e^{-\frac{1}{2}y^{2}}\,dy$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space with a filtration $\mathscr{F}=(\mathscr{F}_{t})_{t\in[0,T]}$ that satisfies the usual conditions, let $W\colon[0,T]\times\Omega\to\mathbb{R}$ be a standard $(\Omega,\mathcal{F},\mathbb{P},\mathscr{F})$-Brownian motion with continuous sample paths, and let $X\colon[0,T]\times\Omega\to\mathbb{R}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

Then it holds for all $K\in\mathbb{R}$ that

[TABLE]

4.1.2 Approximating American options with Bermudan options

In our numerical simulations, we approximate Bermudan options with a finite number of execution times rather than American options, which theoretically can be executed at infinitely many time points (any time before maturity). However, the following result shows that the prices of American options can be approximated with prices of Bermudan options with equidistant execution times if the number of execution times is sufficiently large.

Lemma 4.5**.**

Let $T\in(0,\infty)$, $d,N\in\mathbb{N}$, $\mathfrak{c},\mathfrak{C}\in\mathbb{R}$, let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space with a filtration $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$, let $X\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be a stochastic process with continuous sample paths, let $\left\|\cdot\right\|\colon\mathbb{R}^{d}\to[0,\infty)$ be a norm, assume $\mathbb{E}[\sup_{t\in[0,T]}\|X_{t}\|]<\infty$, and let $\lceil\cdot\rceil\colon[0,T]\to[0,T]$ and $g\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ satisfy for all $t,\mathfrak{t}\in[0,T]$, $x,\mathfrak{x}\in\mathbb{R}^{d}$ that $\lceil t\rceil=\min([t,T]\cap\{0,\nicefrac{{T}}{{N}},\dots,T\})$ and $|g(t,x)-g(\mathfrak{t},\mathfrak{x})|\leq\mathfrak{c}|t-\mathfrak{t}|^{\nicefrac{{1}}{{2}}}+\mathfrak{C}\|x-\mathfrak{x}\|$. Then

1. (i)

it holds that \mathbb{E}\bigl{[}\sup_{t\in[0,T]}|g(t,X_{t})|\bigr{]}\leq|g(0,0)|+\mathfrak{c}T^{\nicefrac{{1}}{{2}}}+\mathfrak{C}\,\mathbb{E}\bigl{[}\sup\nolimits_{t\in[0,T]}\|X_{t}\|\bigr{]}<\infty and 2. (ii)

it holds that

[TABLE]

Proof of

Lemma 4.5.

Throughout this proof let $\lfloor\cdot\rfloor\colon[0,T]\to[0,T]$ satisfy for all $t\in[0,T]$ that $\lfloor t\rfloor=\max([0,t]\cap\{0,\nicefrac{{T}}{{N}},\dots,T\})$. Note that the assumption that for all $t,\mathfrak{t}\in[0,T]$, $x,\mathfrak{x}\in\mathbb{R}^{d}$ it holds that $|g(t,x)-g(\mathfrak{t},\mathfrak{x})|\leq\mathfrak{c}|t-\mathfrak{t}|^{\nicefrac{{1}}{{2}}}+\mathfrak{C}\|x-\mathfrak{x}\|$ implies that

[TABLE]

This establishes (i). Next observe that the fact $\forall\,r\in\{0,\nicefrac{{T}}{{N}},\dots,T\},\,t\in[r,T]\colon r\leq\lfloor t\rfloor$ ensures that for all $r,t\in[0,T]$ with $\lceil r\rceil\leq t$ it holds that $\lceil r\rceil\leq\lfloor t\rfloor$. This implies that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ and every $t\in[0,T]$ it holds that

[TABLE]

The fact that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ and every $t\in[0,T]$ it holds that $\left\{\rho\leq\lfloor t\rfloor\right\}\in\mathbb{F}_{\lfloor t\rfloor}$ and the fact that for every $t\in[0,T]$ it holds that $\mathbb{F}_{\lfloor t\rfloor}\subseteq\mathbb{F}_{t}$ hence show that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ and every $t\in[0,T]$ it holds that $\left\{\lceil\rho\rceil\leq t\right\}\in\mathbb{F}_{t}$. This proves that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ it holds that $\Omega\ni\omega\mapsto\lceil\rho(\omega)\rceil\in[0,T]$ is an $\mathbb{F}$-stopping time. Hence, we obtain that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ it holds that $\Omega\ni\omega\mapsto\lceil\rho(\omega)\rceil\in\{0,\nicefrac{{T}}{{N}},\dots,T\}$ is an $(\mathbb{F}_{\nicefrac{{nT}}{{N}}})_{n\in\{0,1,\dots,N\}}$-stopping time. This ensures that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ it holds that

[TABLE]

Therefore, we obtain that for every $\mathbb{F}$-stopping time $\rho\colon\Omega\to[0,T]$ it holds that

[TABLE]

This implies that

[TABLE]

This shows (ii). The proof of Lemma 4.5 is thus complete. ∎

4.2 Setting

Framework 4.6**.**

Assume Framework 3.2, let $\zeta_{1}=0.9$, $\zeta_{2}=0.999$, $\varepsilon\in(0,\infty)$, $(\gamma_{m})_{m\in\mathbb{N}}\subseteq(0,\infty)$, $\xi=(\xi_{1},\ldots,\xi_{d})\in\mathbb{R}^{d}$, let $\mathscr{F}=(\mathscr{F}_{t})_{t\in[0,T]}$ be a filtration on $(\Omega,\mathcal{F},\mathbb{P})$ that satisfies the usual conditions, let $W^{m,j}=(W^{m,j,(1)},\ldots,W^{m,j,(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$, $m\in\mathbb{N}_{0}$, $j\in\mathbb{N}$, be independent standard $(\Omega,\mathcal{F},\mathbb{P},\mathscr{F})$-Brownian motions with continuous sample paths, let $\mu\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $\sigma\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ be Lipschitz continuous functions, let $X=(X^{(1)},\ldots,X^{(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

assume for all $n\in\{0,1,\ldots,N\}$ that $\varrho=2\nu$, $\Xi_{0}=0$, and $t_{n}=\frac{nT}{N}$, and assume for all $m\in\mathbb{N}$, $x=(x_{1},\ldots,x_{\nu}),\,y=(y_{1},\ldots,y_{\nu}),\,\eta=(\eta_{1},\ldots,\eta_{\nu})\in\mathbb{R}^{\nu}$ that

[TABLE]

and

[TABLE]

Equations (97)–(98) in Framework 4.6 describe the Adam optimiser with possibly varying learning rates (cf. Kingma & Ba [60], e.g., E, Han, & Jentzen [37, (4.3)–(4.4) in Subsection 4.1 and (5.4)–(5.5) in Subsection 5.2], and lines 93–96 in Python code 1 in Subsection 5.1 below). Furthermore, in the context of pricing American-style financial derivatives, we think

• •

of $T$ as the maturity,

• •

of $d$ as the dimension of the associated optimal stopping problem,

• •

of $N$ as the time discretisation parameter employed,

• •

of $M$ as the total number of training steps employed in the Adam optimiser,

• •

of $g$ as the discounted pay-off function,

• •

of $\{t_{0},t_{1},\ldots,t_{N}\}$ as the discrete time grid employed,

• •

of $J_{0}$ as the number of Monte Carlo samples employed in the final integration for the price approximation,

• •

of $(J_{m})_{m\in\mathbb{N}}$ as the sequence of batch sizes employed in the Adam optimiser,

• •

of $\zeta_{1}$ as the momentum decay factor, of $\zeta_{2}$ as the second momentum decay factor, and of $\varepsilon$ as the regularising factor employed in the Adam optimiser,

• •

of $(\gamma_{m})_{m\in\mathbb{N}}$ as the sequence of learning rates employed in the Adam optimiser,

• •

and, where applicable, of $X$ as a continuous-time model for $d$ underlying stock prices with initial prices $\xi$, drift coefficient function $\mu$, and diffusion coefficient function $\sigma$.

Moreover, note that for every $m\in\mathbb{N}_{0}$, $j\in\mathbb{N}$ the stochastic processes $W^{m,j,(1)}=\linebreak(W^{m,j,(1)}_{t})_{t\in[0,T]}$, …, $W^{m,j,(d)}=(W^{m,j,(d)}_{t})_{t\in[0,T]}$ are the components of the $d$-dimensional standard Brownian motion $W^{m,j}=(W^{m,j}_{t})_{t\in[0,T]}$ and hence constitute each a one-dimensional standard Brownian motion.

4.3 Examples with known one-dimensional representation

In this subsection, we test the algorithm of Framework 3.2 on different $d$-dimensional optimal stopping problems that can be represented as one-dimensional optimal stopping problems. This representation allows us to employ a numerical method for the one-dimensional optimal stopping problem to compute reference values for the original $d$-dimensional optimal stopping problem. We refer to Subsection 4.4 below for more challenging examples where a one-dimensional representation is not known.

4.3.1 Optimal stopping of a Brownian motion

4.3.1.1 A Bermudan put-type example with three exercise opportunities

In this subsection, we test the algorithm of Framework 3.2 on the example of optimally stopping a correlated Brownian motion under a put option inspired pay-off function with three possible exercise dates. Among other things, we examine the performance of the algorithm for different numbers of hidden layers of the employed neural networks.

Assume Framework 4.6, let $H\in\mathbb{N}_{0}$, $r=0.02=2\%$, $\beta=0.3=30\%$, $\chi=95$, $K=90$, $Q=(Q_{i,j})_{(i,j)\in\{1,\ldots,d\}^{2}},\mathfrak{S}\in\mathbb{R}^{d\times d}$ satisfy for all $i\in\{1,\ldots,d\}$ that $Q_{i,i}=1$, $\forall\,j\in\{1,\ldots,d\}\setminus\{i\}\colon Q_{i,j}=0.1$, and $\mathfrak{S}\,\mathfrak{S}^{*}=Q$, let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $W^{0,1}$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $W^{0,1,(1)}$, assume that each of the employed neural networks has $H$ hidden layers, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,2\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=1$, $N=2$, $M=500$, $\mathcal{X}_{n}^{m-1,j}=\mathfrak{S}\,W_{t_{n}}^{m-1,j}$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,100]}(m)+10^{-3}\,\mathbbm{1}_{(100,300]}(m)+10^{-4}\,\mathbbm{1}_{(300,\infty)}(m)]$, and

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the real number

[TABLE]

The numbers in Table 1 were obtained with our standard network architecture with two hidden layers. It shows approximations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$, and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $d\in\{1,5,10,50,100,500,1000\}$. For each case, the calculations of the results in Tables 1–3 are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python. Furthermore, in the approximative calculations of the relative approximation error associated with $\mathcal{P}$, the exact number (100) was replaced, independently of the dimension $d$, by the real number

[TABLE]

(cf. Corollary 4.2), which, in turn, was replaced by the value $7.894$. The latter was computed in Matlab R2017b using the binomial tree method implemented as Matlab’s function optstockbycrr with 20,000 nodes. Note that (101) corresponds to the price of a Bermudan put option on a single stock in the Black–Scholes model with initial stock price $\chi$, interest rate $r$, volatility $\beta$, strike price $K$, maturity $T$, and $N$ possible exercise dates.

Due to the underlying one-dimensional structure, all examples in Subsection 4.3 admit an optimal stopping rule which, at each possible exercise date $t_{n}$, checks whether the current pay-off is above a threshold $c_{n}\in\mathbb{R}$. Therefore, we also apply our algorithm to the example in Subsection 4.3.1.1 using networks with one input neuron and no hidden layers. This corresponds to learning the thresholds $c_{n}$ from simulated pay-offs with one-dimensional logistic regressions. We used the same number of simulations as in Table 1 and batch normalisation before the first linear transformation but no batch normalisation before the logistic function. As can be seen from Table 2, the resulting approximations have the same accuracy as the ones of Table 1 and, in addition, the computations times are shorter. However, as we will see in Subsection 4.4, it cannot be hoped that good results can be obtained with a simplified network architecture if the stopping problem is more complex.

Table 3 shows approximations of (100) for $d=10$ obtained with networks with $(d+1)$-dimensional input layers and different numbers of hidden layers. Again, it can be seen that in this example hidden layers do not improve the accuracy of the results.

4.3.1.2 An American put-type example

In this subsection, we test the algorithm of Framework 3.2 on the example of optimally stopping a standard Brownian motion under a put option inspired pay-off function.

Assume Framework 4.6, let $r=0.06=6\%$, $\beta=0.4=40\%$, $\chi=K=40$, let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $W^{0,1}$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $W^{0,1,(1)}$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=1$, $N=50$, $M=1500\,\mathbbm{1}_{[1,50]}(d)+1800\,\mathbbm{1}_{(50,100]}(d)+3000\,\mathbbm{1}_{(100,\infty)}(d)$, $\mathcal{X}_{n}^{m-1,j}=W_{t_{n}}^{m-1,j}$, $J_{0}=$ 4,096,000, $J_{m}=8192\,\mathbbm{1}_{[1,50]}(d)+4096\,\mathbbm{1}_{(50,100]}(d)+2048\,\mathbbm{1}_{(100,\infty)}(d)$, $\varepsilon=0.001$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,\nicefrac{{M}}{{3}}]}(m)+10^{-3}\,\mathbbm{1}_{(\nicefrac{{M}}{{3}},\nicefrac{{2M}}{{3}}]}(m)+10^{-4}\,\mathbbm{1}_{(\nicefrac{{2M}}{{3}},\infty)}(m)]$, and

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the real number

[TABLE]

We report approximations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$, and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $d\in\{1,5,10,50,100,500,1000\}$ in Table 4. For each case, the calculations of the results in Table 4 are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python. Furthermore, in the approximative calculations of the relative approximation error associated with $\mathcal{P}$, the exact number (103) was replaced, independently of the dimension $d$, by the real number

[TABLE]

(cf. Corollary 4.2), which, in turn, was replaced by the value $5.318$ (cf. Longstaff & Schwartz [75, Table 1 in Section 3]). This value was calculated using the binomial tree method on Smirnov’s website [87] with 20,000 nodes. Note that (104) corresponds to the price of an American put option on a single stock in the Black–Scholes model with initial stock price $\chi$, interest rate $r$, volatility $\beta$, strike price $K$, and maturity $T$. Moreover, note that in this example we have that for all $m\in\mathbb{N}$ the batch size is $J_{m}=8192\,\mathbbm{1}_{[1,50]}(d)+4096\,\mathbbm{1}_{(50,100]}(d)+2048\,\mathbbm{1}_{(100,\infty)}(d)$. In particular, for all $m\in\mathbb{N}$ the batch size $J_{m}$ is decreasing as $d$ increases. Such kind of fine-tuning is necessary to fit the examples into the memory available on the GPU here and below.

4.3.2 Geometric average-type options

4.3.2.1 An American geometric average put-type example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing an American geometric average put-type option on up to 200 distinguishable stocks in the Black–Scholes model. Moreover, we compare its performance with that of the algorithm introduced in [9].

Assume Framework 4.6, assume that $d\in\{40,80,120,\ldots\}$, let $\beta=(\beta_{1},\ldots,\beta_{d})\in\mathbb{R}^{d}$, $\rho,\tilde{\delta},\tilde{\beta},\delta_{1},\delta_{2},\ldots,\delta_{d}\in\mathbb{R}$, $r=0.6$, $K=95$, $\tilde{\xi}=100$ satisfy for all $i\in\{1,\ldots,d\}$ that $\beta_{i}=\min\{0.04\,[(i-1)\bmod 40],1.6-0.04\,[(i-1)\bmod 40]\}$, $\rho=\frac{1}{d}\|\beta\|_{\mathbb{R}^{d}}^{2}=\frac{1}{40}\sum_{i=1}^{40}(\beta_{i})^{2}=0.2136$, \delta_{i}=r-\frac{\rho}{d}\bigl{(}i-\frac{1}{2}\bigr{)}-\frac{1}{5\sqrt{d}}, $\tilde{\delta}=r-\frac{1}{\sqrt{d}}\sum_{i=1}^{d}(r-\delta_{i})+\frac{\sqrt{d}-1}{2d}\|\beta\|_{\mathbb{R}^{d}}^{2}=r-\frac{\rho}{2}-\frac{1}{5}=0.2932$, and $\tilde{\beta}=\tfrac{1}{\sqrt{d}}\|\beta\|_{\mathbb{R}^{d}}=\sqrt{\rho}$, let $Y\colon[0,T]\times\Omega\to\mathbb{R}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $Y$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=1$, $N=100$, $M=1800\,\mathbbm{1}_{[1,120]}(d)+3000\,\mathbbm{1}_{(120,\infty)}(d)$, $J_{0}=$ 4,096,000, $J_{m}=8192\,\mathbbm{1}_{[1,120]}(d)+4096\,\mathbbm{1}_{(120,\infty)}(d)$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,\nicefrac{{M}}{{3}}]}(m)+10^{-3}\,\mathbbm{1}_{(\nicefrac{{M}}{{3}},\nicefrac{{2M}}{{3}}]}(m)+10^{-4}\,\mathbbm{1}_{(\nicefrac{{2M}}{{3}},\infty)}(m)]$, $\xi_{i}=(100)^{\nicefrac{{1}}{{\sqrt{d}}}}$, $\mu(x)=((r-\delta_{1})\,x_{1},\ldots,(r-\delta_{d})\,x_{d})$, $\sigma(x)=\operatorname{diag}(\beta_{1}x_{1},\ldots,\beta_{d}x_{d})$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

In Table 5, we show approximations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$, and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $d\in\{40,80,120,160,200\}$. For each case, the calculations of the results in Table 5 are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python. Furthermore, in the approximative calculations of the relative approximation error associated with $\mathcal{P}$, the exact value of the price (108) was replaced, independently of the dimension $d$, by the real number

[TABLE]

(cf. Proposition 4.3), which, in turn, was replaced by the value $6.545$. The latter was calculated using the binomial tree method on Smirnov’s website [87] with 20,000 nodes. Note that (109) corresponds to the price of an American put option on a single stock in the Black–Scholes model with initial stock price $\tilde{\xi}$, interest rate $r$, dividend yield $\tilde{\delta}$, volatility $\tilde{\beta}$, strike price $K$, and maturity $T$.

In addition, we computed approximations of the price (108) using the algorithm introduced in [9], where the random variable $\hat{L}$ from [9, Subsection 3.1] plays the role analogous to $\mathcal{P}$. In order to maximise comparability, the hyperparameters and neural network architectures employed for the algorithm from [9] were chosen to be identical to the corresponding ones used for computing realisations of $\mathcal{P}$. Table 6 shows approximations of the mean of $\hat{L}$ (cf. [9, Subsection 3.1]), of the standard deviation of $\hat{L}$, and of the relative $L^{1}$-approximation error associated with $\hat{L}$, the uncorrected sample standard deviation of the relative approximation error associated with $\hat{L}$, and the average runtime in seconds needed for calculating one realisation of $\hat{L}$ for $d\in\{40,80,120,160,200\}$. For each case, the calculations of the results in Table 6 are based on $10$ independent realisations of $\hat{L}$, which were obtained from an implementation in Python. In the approximative calculations of the relative approximation error associated with $\hat{L}$, the exact value of the price (108) again was replaced, independently of the dimension $d$, by the value $6.545$. Comparing Table 5 with Table 6, we note that in the present cases the algorithm of Framework 3.2 and the algorithm from [9] exhibit very similar performance in terms of both accuracy and speed, with a slight runtime advantage for the algorithm of Framework 3.2.

4.3.2.2 An American geometric average call-type example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing an American geometric average call-type option on up to 200 correlated stocks in the Black–Scholes model. This example is taken from Sirignano & Spiliopoulos [86, Subsection 4.3].

Assume Framework 4.6, let $r=0\%$, $\delta=0.02=2\%$, $\beta=0.25=25\%$, $K=\tilde{\xi}=1$, $Q=(Q_{i,j})_{(i,j)\in\{1,\ldots,d\}^{2}},\mathfrak{S}=(\varsigma_{1},\ldots,\varsigma_{d})\in\mathbb{R}^{d\times d}$, $\tilde{\delta},\tilde{\beta}\in\mathbb{R}$ satisfy for all $i\in\{1,\ldots,d\}$ that $Q_{i,i}=1$, $\forall\,j\in\{1,\ldots,d\}\setminus\{i\}\colon Q_{i,j}=0.75$, $\mathfrak{S}^{*}\,\mathfrak{S}=Q$, $\tilde{\delta}=\delta+\frac{1}{2}(\beta^{2}-(\tilde{\beta})^{2})$, and $\tilde{\beta}=\tfrac{\beta}{2d}\sqrt{d\,(3d+1)}$, let $Y\colon[0,T]\times\Omega\to\mathbb{R}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $Y$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=2$, $N=50$, $M=1600$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,400]}(m)+10^{-3}\,\mathbbm{1}_{(400,800]}(m)\allowbreak+10^{-4}\,\mathbbm{1}_{(800,\infty)}(m)]$, $\xi_{i}=1$, $\mu(x)=(r-\delta)\,x$, $\sigma(x)=\beta\operatorname{diag}(x_{1},\dots,x_{d})\,\mathfrak{S}^{*}$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

We report approximations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, of the real number

[TABLE]

and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$, and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $d\in\{3,20,100,200\}$ in Table 7. The approximative calculations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the computations of the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$ as well as the computations of the average runtime for calculating one realisation of $\mathcal{P}$ in Table 7 each are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python. Furthermore, in the approximative calculations of the relative approximation error associated with $\mathcal{P}$, the exact value of the price (113) was replaced by the number (114) (cf. Proposition 4.3), which was approximatively calculated using the binomial tree method on Smirnov’s website [87] with 20,000 nodes. Note that (114) corresponds to the price of an American call option on a single stock in the Black–Scholes model with initial stock price $\tilde{\xi}$, interest rate $r$, dividend yield $\tilde{\delta}$, volatility $\tilde{\beta}$, strike price $K$, and maturity $T$.

4.3.2.3 Another American geometric average call-type example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing an American geometric average call-type option on up to 400 distinguishable stocks in the Black–Scholes model.

Assume Framework 4.6, assume that $d\in\{40,80,120,\ldots\}$, let $\beta=(\beta_{1},\ldots,\beta_{d})\in\mathbb{R}^{d}$, $\alpha_{1},\ldots,\alpha_{d}\in\mathbb{R}$, $r,\tilde{\beta}\in(0,\infty)$, $K=95$, $\tilde{\xi}=100$ satisfy for all $i\in\{1,\ldots,d\}$ that $\beta_{i}=\frac{0.4\,i}{d}$, $\alpha_{i}=\min\{0.01\,[(i-1)\bmod 40],0.4-0.01\,[(i-1)\bmod 40]\}$, r=\frac{1}{d}\sum_{i=1}^{d}\alpha_{i}-\frac{d-1}{2d^{2}}\|\beta\|_{\mathbb{R}^{d}}^{2}=0.1-\frac{0.08}{d^{2}}(d-1)\bigl{(}\frac{d}{3}+\frac{1}{2}+\frac{1}{6d}\bigr{)}, and \tilde{\beta}=\tfrac{1}{d}\|\beta\|_{\mathbb{R}^{d}}=\frac{0.4}{d}\bigl{(}\frac{d}{3}+\frac{1}{2}+\frac{1}{6d}\bigr{)}^{\nicefrac{{1}}{{2}}}, let $Y\colon[0,T]\times\Omega\to\mathbb{R}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $Y$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=3$, $N=50$, $M=1500$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,\nicefrac{{M}}{{3}}]}(m)+10^{-3}\,\mathbbm{1}_{(\nicefrac{{M}}{{3}},\nicefrac{{2M}}{{3}}]}(m)+10^{-4}\,\mathbbm{1}_{(\nicefrac{{2M}}{{3}},\infty)}(m)]$, $\xi_{i}=100$, $\mu(x)=(\alpha_{1}x_{1},\ldots,\alpha_{d}x_{d})$, $\sigma(x)=\operatorname{diag}(\beta_{1}x_{1},\ldots,\beta_{d}x_{d})$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

In Table 8, we show approximations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, of the real number

[TABLE]

and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$, and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $d\in\{40,80,120,160,200,\allowbreak 400\}$. The approximative calculations of the mean of $\mathcal{P}$, of the standard deviation of $\mathcal{P}$, and of the relative $L^{1}$-approximation error associated with $\mathcal{P}$, the computations of the uncorrected sample standard deviation of the relative approximation error associated with $\mathcal{P}$ as well as the computations of the average runtime for calculating one realisation of $\mathcal{P}$ in Table 8 each are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python. Moreover, in the approximative calculations of the relative approximation error associated with $\mathcal{P}$, the exact value of the price (118) was replaced by the real number

[TABLE]

(cf. Proposition 4.3). It is well known (cf., e.g., Shreve [84, Corollary 8.5.3]) that the number (120) is equal to the number (119), which was approximatively computed in Matlab R2017b using Lemma 4.4 above. Note that (120) corresponds to the price of an American call option on a single stock in the Black–Scholes model with initial stock price $\tilde{\xi}$, interest rate $r$, volatility $\tilde{\beta}$, strike price $K$, and maturity $T$, while (119) corresponds to the price of a European call option on a single stock in the Black–Scholes model with initial stock price $\tilde{\xi}$, interest rate $r$, volatility $\tilde{\beta}$, strike price $K$, and maturity $T$.

4.4 Examples without known one-dimensional representation

In Subsection 4.3 above, numerical results for examples with a one-dimensional representation can be found. We test in this subsection different examples where such a representation is not known.

4.4.1 Max-call options

4.4.1.1 A Bermudan max-call standard benchmark example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing a Bermudan max-call option on up to 500 stocks in the Black–Scholes model (cf. [9, Subsection 4.1]). In the case of up to five underlying stocks, this example is a standard benchmark example in the literature (cf., e.g., [26, Subsection 5.4], [75, Subsection 8.1], [3, Section 4], [51, Subsection 5.1], [80, Subsection 4.3], [41, Subsection 3.9], [19, Subsection 4.2], [24, Subsection 5.3], [15, Subsection 6.1], [54, Subsection 4.1], [82, Subsection 7.2], [17, Subsection 6.1], [73, Subsection 5.2.1]).

Assume Framework 4.6, let $H\in\mathbb{N}_{0}$, $r=0.05=5\%$, $\delta=0.1=10\%$, $\beta=0.2=20\%$, $K=100$, let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, assume that each of the employed neural networks has $H$ hidden layers, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=3$, $N=9$, $M=3000+d$, $J_{0}=$ 4,096,000 $J_{m}=8192$, $\varepsilon=0.1$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,500+\nicefrac{{d}}{{5}}]}(m)+10^{-3}\,\mathbbm{1}_{(500+\nicefrac{{d}}{{5}},1500+\nicefrac{{3d}}{{5}}]}(m)+10^{-4}\,\mathbbm{1}_{(1500+\nicefrac{{3d}}{{5}},\infty)}(m)]$, $\xi_{i}=\xi_{1}$, $\mu(x)=(r-\delta)\,x$, $\sigma(x)=\beta\operatorname{diag}(x_{1},\dots,x_{d})$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

In Table 9, we show approximations of the mean and of the standard deviation of $\mathcal{P}$, binomial approximations as well as 95% confidence intervals for the price (123) according to Andersen & Broadie [3, Table 2 in Section 4] (where available), 95% confidence intervals for the price (123) according to Broadie & Cao [24, Table 3 in Subsection 5.3] (where available), and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $(d,\xi_{1})\in\{2,3,5\}\times\{90,100,110\}$ and $H=2$. The approximative calculations of the mean and of the standard deviation of $\mathcal{P}$ as well as the computations of the average runtime for calculating one realisation of $\mathcal{P}$ in Tables 9–11 each are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python (cf. Python code 2 in Subsection 5.2 below).

In Table 10, we list approximations of the mean and of the standard deviation of $\mathcal{P}$ and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $(d,\xi_{1})\in\{10,20,30,50,100,200,500\}\times\{90,100,110\}$ and $H=2$.

To see the impact of the number of hidden layers used in the neural networks, we additionally report in Table 11 approximation results for $d=5$, $\xi_{1}=100$, and $H\in\{0,1,2,3,4,5\}$. We used the same number of simulations as in Tables  9–10, but due to the higher number of hidden layers, we chose $M=5000$ and $\forall\,m\in\mathbb{N}\colon\gamma_{m}=5[10^{-2}\mathbbm{1}_{[1,1000]}(m)+10^{-3}\mathbbm{1}_{(1000,3000]}(m)+10^{-4}\mathbbm{1}_{(3000,\infty)}(m)]$. It can be seen that, in this example, two hidden layers yield better results than zero or one hidden layer. But more than two hidden layers do not lead to an improvement.

4.4.1.2 A high-dimensional Bermudan max-call benchmark example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing the Bermudan max-call option from the example in Subsection 4.4.1.1 in a case with 5000 underlying stocks. All Python source codes corresponding to this example were run in single precision (float32) on a NVIDIA Tesla P100 GPU.

Assume Framework 4.6, let $r=0.05=5\%$, $\delta=0.1=10\%$, $\beta=0.2=20\%$, $K=100$, let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=3$, $d=5000$, $N=9$, $J_{0}=2^{20}$, $J_{m}=1024$, $\varepsilon=10^{-8}$, $\gamma_{m}=10^{-2}\,\mathbbm{1}_{[1,2000]}(m)+10^{-3}\,\mathbbm{1}_{(2000,4000]}(m)+10^{-4}\,\mathbbm{1}_{(4000,\infty)}(m)$, $\xi_{i}=100$, $\mu(x)=(r-\delta)\,x$, $\sigma(x)=\beta\operatorname{diag}(x_{1},\dots,x_{d})$, that

[TABLE]

and that

[TABLE]

For sufficiently large $M\in\mathbb{N}$, the random variable $\mathcal{P}$ provides approximations of the price

[TABLE]

In Table 12, we show a realisation of $\mathcal{P}$, a 95% confidence interval for the corresponding realisation of the random variable

[TABLE]

the corresponding realisation of the relative approximation error associated with $\mathcal{P}$, and the runtime in seconds needed for calculating the realisation of $\mathcal{P}$ for $M\in\{0,250,500,\ldots,\allowbreak 2000\}\cup\{6000\}$. In addition, Figure 1 depicts a realisation of the relative approximation error associated with $\mathcal{P}$ against $M\in\{0,10,20,\ldots,2000\}$. For each case, the 95% confidence interval for the realisation of the random variable (127) in Table 12 was computed based on the corresponding realisation of $\mathcal{P}$, the corresponding sample standard deviation, and the $0.975$ quantile of the standard normal distribution (cf., e.g., [9, Subsection 3.3]). Moreover, in the approximative calculations of the realisation of the relative approximation error associated with $\mathcal{P}$, in Table 12 and Figure 1 the exact value of the price (126) was replaced by the value 165.430, which corresponds to a realisation of $\mathcal{P}$ with $M=6000$ (cf. Table 12). All realisations of $\mathcal{P}$ have been obtained from running a close variant of the Python code used for the example in Subsection 4.4.1.1 with $d=5000$ (cf. Python code 2 in Subsection 5.2 below).

4.4.1.3 Another Bermudan max-call example

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing a Bermudan max-call option for different maturities and strike prices on up to 400 correlated stocks, that do not pay dividends, in the Black–Scholes model. This example is taken from Barraquand & Martineau [5, Section VII].

Assume Framework 4.6, let $\eta=\nicefrac{{30}}{{365}}$, $r=0.05\,\eta=5\%\cdot\eta$, $\beta=0.4\sqrt{\eta}=40\%\cdot\sqrt{\eta}$, $K\in\{35,40,45\}$, $Q=(Q_{i,j})_{(i,j)\in\{1,\ldots,d\}^{2}},\mathfrak{S}=(\varsigma_{1},\ldots,\varsigma_{d})\in\mathbb{R}^{d\times d}$ satisfy for all $i\in\{1,\ldots,d\}$ that $Q_{i,i}=1$, $\forall\,j\in\{1,\ldots,d\}\setminus\{i\}\colon Q_{i,j}=0.5$, and $\mathfrak{S}^{*}\,\mathfrak{S}=Q$, let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $N=10$, $M=1600$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=0.001$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,400]}(m)+10^{-3}\,\mathbbm{1}_{(400,800]}(m)+10^{-4}\,\mathbbm{1}_{(800,\infty)}(m)]$, $\xi_{i}=40$, $\mu(x)=r\,x$, $\sigma(x)=\beta\operatorname{diag}(x_{1},\dots,x_{d})\,\mathfrak{S}^{*}$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

In Table 13, we show approximations of the mean and of the standard deviation of $\mathcal{P}$, Monte Carlo approximations of the European max-call option price

[TABLE]

corresponding to (130), approximations of the price (130) according to [5, Table 4 in Section VII] (where available), and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for

[TABLE]

The approximative calculations of the mean and of the standard deviation of $\mathcal{P}$ as well as the computations of the average runtime for calculating one realisation of $\mathcal{P}$ in Table 13 each are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python (cf. Python code 3 in Subsection 5.3 below). Furthermore, the Monte Carlo approximations of the European price (131) in Table 13 each were calculated in double precision (float64) and are based on $2\cdot 10^{10}$ independent realisations of the random variable $\Omega\ni\omega\mapsto g(T,X_{T}(\omega))\in\mathbb{R}$.

4.4.2 A strangle spread basket option

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing an American strangle spread basket option on five correlated stocks in the Black–Scholes model. This example is taken from Kohler, Krzyżak, & Todorovic [66, Section 4] (cf. also Kohler [63, Section 3] and Kohler, Krzyżak, & Walk [67, Section 4]).

Assume Framework 4.6, let $r=0.05=5\%$, $K_{1}=75$, $K_{2}=90$, $K_{3}=110$, $K_{4}=125$, let $\mathfrak{S}=(\varsigma_{1},\ldots,\varsigma_{5})\in\mathbb{R}^{5\times 5}$ be given by

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $X$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $i\in\{1,\ldots,d\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=1$, $d=5$, $N=48$, $M=750$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,250]}(m)+10^{-3}\,\mathbbm{1}_{(250,500]}(m)+10^{-4}\,\mathbbm{1}_{(500,\infty)}(m)]$, $\xi_{i}=100$, $\mu(x)=r\,x$, $\sigma(x)=\operatorname{diag}(x_{1},\dots,x_{d})\,\mathfrak{S}^{*}$, that

[TABLE]

and that

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

Table 14 shows approximations of the mean and of the standard deviation of $\mathcal{P}$, a lower bound for the price (136) according to Kohler, Krzyżak, & Todorovic [66, Figure 4.5 in Section 4] (cf. also Kohler [63, Figure 2 in Section 3] and, for an upper bound for the price (136), Kohler, Krzyżak, & Walk [67, Figure 4.2 in Section 4]), and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$. Since the mean of $\mathcal{P}$ is also a lower bound for the price (136), a higher value indicates a better approximation of the price (136) (cf. Table 14). The approximative calculations of the mean and of the standard deviation of $\mathcal{P}$ as well as the computation of the average runtime for calculating one realisation of $\mathcal{P}$ in Table 14 each are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python (cf. Python code 4 in Subsection 5.4 below).

4.4.3 A put basket option in Dupire’s local volatility model

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing an American put basket option on five stocks in Dupire’s local volatility model. This example is taken from Labart & Lelong [70, Subsection 6.3] with the modification that we also consider the case where the underlying stocks do not pay any dividends.

Assume Framework 4.6, let $L=10$, $r=0.05=5\%$, $\delta\in\{0\%,10\%\}$, $K=100$, assume for all $i\in\{1,\ldots,d\}$, $x\in\mathbb{R}^{d}$ that $\xi_{i}=100$ and $\mu(x)=(r-\delta)\,x$, let $\beta\colon[0,T]\times\mathbb{R}\to\mathbb{R}$ and $\boldsymbol{\sigma}\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ be the functions which satisfy for all $t\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that

[TABLE]

and $\boldsymbol{\sigma}(t,x)=\operatorname{diag}(\beta(t,x_{1}),\beta(t,x_{2}),\dots,\beta(t,x_{d}))$, let $S=(S^{(1)},\ldots,S^{(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$ be an $\mathscr{F}$-adapted stochastic process with continuous sample paths which satisfies that for all $t\in[0,T]$ it holds $\mathbb{P}$-a.s. that

[TABLE]

let $\mathcal{Y}^{m,j}=(\mathcal{Y}^{m,j,(1)},\ldots,\mathcal{Y}^{m,j,(d)})\colon[0,T]\times\Omega\to\mathbb{R}^{d}$, $j\in\mathbb{N}$, $m\in\mathbb{N}_{0}$, be the stochastic processes which satisfy for all $m\in\mathbb{N}_{0}$, $j\in\mathbb{N}$, $\ell\in\{0,1,\dots,L-1\}$, t\in\bigl{[}\frac{\ell T}{L},\frac{(\ell+1)T}{L}\bigr{]}, $i\in\{1,\ldots,d\}$ that $\mathcal{Y}_{0}^{m,j,(i)}=\log(\xi_{i})$ and

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{t})_{t\in[0,T]}$ be the filtration generated by $S$, let $\mathfrak{F}=(\mathfrak{F}_{t})_{t\in[0,T]}$ be the filtration generated by $\mathcal{Y}^{0,1}$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T=1$, $d=5$, $M=1200$, $\mathcal{X}_{n}^{m-1,j}=\mathcal{Y}_{t_{n}}^{m-1,j}$, $J_{0}=$ 4,096,000, $J_{m}=8192$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,400]}(m)+10^{-3}\,\mathbbm{1}_{(400,800]}(m)+10^{-4}\,\mathbbm{1}_{(800,\infty)}(m)]$, and

[TABLE]

The random variable $\mathcal{P}$ given in (67) provides approximations of the price

[TABLE]

which, in turn, is an approximation of the price

[TABLE]

In Table 15, we show approximations of the mean and of the standard deviation of $\mathcal{P}$ and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $(\delta,N)\in\{0\%,10\%\}\times\{5,10,50,100\}$. For each case, the calculations of the results in Table 15 are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python (cf. Python code 5 in Subsection 5.5 below). According to [70, Subsection 6.3], the value $6.30$ is an approximation of the price (141) for $\delta=10\%$. Furthermore, the European put basket option price \mathbb{E}\bigl{[}g(T,\mathcal{Y}_{T}^{0,1})\bigr{]} corresponding to (141) was approximatively calculated using a Monte Carlo approximation based on $10^{10}$ realisations of the random variable $\Omega\ni\omega\mapsto g(T,\mathcal{Y}_{T}^{0,1}(\omega))\in\mathbb{R}$ (cf. Python code 5 in Subsection 5.5 below), which resulted in the value $1.741$ in the case $\delta=0\%$ and in the value $6.304$ in the case $\delta=10\%$.

4.4.4 A path-dependent financial derivative

In this subsection, we test the algorithm of Framework 3.2 on the example of pricing a specific path-dependent financial derivative contingent on prices of a single underlying stock in the Black–Scholes model, which is formulated as a 100-dimensional optimal stopping problem. This example is taken from Tsitsiklis & Van Roy [90, Section IV] with the modification that we consider a finite instead of an infinite time horizon.

Assume Framework 4.6, let $r=0.0004=0.04\%$, $\beta=0.02=2\%$, let $\mathcal{W}^{m,j}\colon[0,\infty)\times\Omega\to\mathbb{R}$, $j\in\mathbb{N}$, $m\in\mathbb{N}_{0}$, be independent $\mathbb{P}$-standard Brownian motions with continuous sample paths, let $S^{m,j}\colon[-100,\infty)\times\Omega\to\mathbb{R}$, $j\in\mathbb{N}$, $m\in\mathbb{N}_{0}$, and $\mathcal{Y}^{m,j}\colon\mathbb{N}_{0}\times\Omega\to\mathbb{R}^{100}$, $j\in\mathbb{N}$, $m\in\mathbb{N}_{0}$, be the stochastic processes which satisfy for all $m,n\in\mathbb{N}_{0}$, $j\in\mathbb{N}$, $t\in[-100,\infty)$ that S_{t}^{m,j}=\exp\bigl{(}\bigl{[}r-\tfrac{1}{2}\beta^{2}\bigr{]}(t+100)+\beta\,\mathcal{W}_{t+100}^{m,j}\bigr{)}\,\xi_{1} and

[TABLE]

let $\mathbb{F}=(\mathbb{F}_{n})_{n\in\mathbb{N}_{0}}$ be the filtration generated by $\mathcal{Y}^{0,1}$, and assume for all $m,j\in\mathbb{N}$, $n\in\{0,1,\ldots,N\}$, $s\in[0,T]$, $x=(x_{1},\ldots,x_{d})\in\mathbb{R}^{d}$ that $T\in\mathbb{N}$, $d=100$, $N=T$, $M=1200\,\mathbbm{1}_{[1,150]}(T)+1500\,\mathbbm{1}_{(150,250]}(T)+3000\,\mathbbm{1}_{(250,\infty)}(T)$, $\mathcal{X}_{n}^{m-1,j}=\mathcal{Y}_{n}^{m-1,j}$, $J_{0}=$ 4,096,000, $J_{m}=8192\,\mathbbm{1}_{[1,150]}(T)+4096\,\mathbbm{1}_{(150,250]}(T)+512\,\mathbbm{1}_{(250,\infty)}(T)$, $\varepsilon=10^{-8}$, $\gamma_{m}=5\,[10^{-2}\,\mathbbm{1}_{[1,\nicefrac{{M}}{{3}}]}(m)+10^{-3}\,\mathbbm{1}_{(\nicefrac{{M}}{{3}},\nicefrac{{2M}}{{3}}]}(m)+10^{-4}\,\mathbbm{1}_{(\nicefrac{{2M}}{{3}},\infty)}(m)]$, and $g(s,x)=e^{-rs}\,x_{100}$. The random variable $\mathcal{P}$ provides approximations of the real number

[TABLE]

In Table 16, we show approximations of the mean and of the standard deviation of $\mathcal{P}$ and the average runtime in seconds needed for calculating one realisation of $\mathcal{P}$ for $T\in\{100,150,200,250,1000\}$. For each case, the calculations of the results in Table 16 are based on $10$ independent realisations of $\mathcal{P}$, which were obtained from an implementation in Python (cf. Python code 6 in Subsection 5.6 below). Note that in this example time is measured in days and that, roughly speaking, (144) corresponds to the price of a financial derivative which, if the holder decides to exercise, pays off the amount given by the ratio between the current underlying stock price and the underlying stock price 100 days ago (cf. [90, Section IV] for more details). According to [90, Subsection IV.D], the value $1.282$ is a lower bound for the price

[TABLE]

which corresponds to the price (144) in the case of an infinite time horizon. Since the mean of $\mathcal{P}$ is a lower bound for the price (144), which, in turn, is a lower bound for the price (145), a higher value indicates a better approximation of the price (145). In addition, observe that the price (144) is non-decreasing in $T$. While in our numerical simulations the approximate value of the mean of $\mathcal{P}$ is less or equal than $1.282$ for comparatively small time horizons, i.e., for $T\leq 150$, it is already higher for slightly larger time horizons, i.e., for $T\geq 200$ (cf. Table 16).

5 Python source codes

In Subsections 5.2–5.6 below we present Python source codes associated to the numerical simulations in Subsections 4.4.1.1, 4.4.1.3, 4.4.2, 4.4.3, and 4.4.4 above. The following Python source code, Python code 1 in Subsection 5.1 below, contains the main part of the implementation of the algorithm in Framework 3.2 and Framework 4.6 and is employed in each of the Python source codes in Subsections 5.2–5.6.

5.1 A Python source code for the algorithm

5.2 A Python source code associated to the numerical simulations in Subsection 4.4.1.1

5.3 A Python source code associated to the numerical simulations in Subsection 4.4.1.3

5.4 A Python source code associated to the numerical simulations in Subsection 4.4.2

5.5 A Python source code associated to the numerical simulations in Subsection 4.4.3

5.6 A Python source code associated to the numerical simulations in Subsection 4.4.4

Acknowledgements

This project has been partially supported through the ETH Research Grant ETH-47 15-2 ‘Mild stochastic calculus and numerical approximations for nonlinear stochastic evolution equations with Lévy noise’, by the project ‘Construction of New Smoothness Spaces on Domains’ (project number I 3403) funded by the Austrian Science Fund (FWF), and by the project ‘Deep artificial neural network approximations for stochastic partial differential equations: Algorithms and convergence proofs’ (project number 184220) funded by the Swiss National Science Foundation (SNSF).