1. Orthogonal transforms and Brownian paths {#s000005}
===========================================

There are several constructions that are frequently used to construct *discrete Brownian paths*, by which we mean a random function $B$ on a given set $\left\{ t_{1},\ldots,t_{n} \right\} \subseteq \mathbb{R}$, $0 < t_{1} < \cdots < t_{n} \leq 1$, so that $B = \left( B_{t_{1}},\ldots,B_{t_{n}} \right)$ is a Gaussian vector with mean zero and covariance matrix $$\left( \min\left( t_{j},t_{k} \right) \right)_{j,k = 1}^{n} = \begin{pmatrix}
t_{1} & t_{1} & t_{1} & \ldots & t_{1} \\
t_{1} & t_{2} & t_{2} & \ldots & t_{2} \\
t_{1} & t_{2} & t_{3} & \ldots & t_{3} \\
 \vdots & \vdots & \vdots & \ddots & \vdots \\
t_{1} & t_{2} & t_{3} & \ldots & t_{n} \\
\end{pmatrix}\text{.}$$

The case where the $t_{j}$ are evenly spaced is the most important one from the practical point of view. In that case the covariance matrix equals $$\left( \frac{1}{n}\min\left( j,k \right) \right)_{j,k = 1}^{n} = \frac{1}{n}\begin{pmatrix}
1 & 1 & 1 & \ldots & 1 \\
1 & 2 & 2 & \ldots & 2 \\
1 & 2 & 3 & \ldots & 3 \\
 \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & 2 & 3 & \ldots & n \\
\end{pmatrix}\text{.}$$Throughout the paper this matrix will be denoted by $\Sigma^{(n)}$ or, if there is no danger of confusion, simply by $\Sigma$.

The arguably most straightforward method (a.k.a. *forward method* or *step-by-step method*) is to compute the cumulative sum of $n$ independent normal variables of mean zero and variance $\frac{1}{n}$. All constructions we present in this article are equivalent to this simple method from the probabilistic point of view.

However, there are refined simulation methods, for example stratified sampling (cf. [@br000050]) and quasi-Monte Carlo methods (see [@br000095]), which achieve higher convergence rates for some problems. Those techniques have in common that they require the identification of more important and less important input variables. For many problems the straightforward method does not provide this.

For this reason alternatives to the forward construction are frequently used, the Brownian bridge (BB) construction (a.k.a. Lévy--Ciesielski construction or midpoint displacement) and the principal component analysis (PCA) construction (a.k.a. singular value construction). The first use of the BB construction in finance is due to [@br000090], the first use of PCA for financial applications was presented in [@br000010], both with dramatic improvement of convergence rates.

It has been mentioned by Papageorgiou [@br000100] that in fact any decomposition $AA^{T} = \Sigma$ provides a construction for a discrete approximation of a Brownian path via $Y = AX$, where $X$ is a standard normal vector. In that context, the forward construction corresponds to the Cholesky decomposition of $\Sigma$, $\Sigma = SS^{\top}$, where $S$ is the summation operator $$S = \frac{1}{\sqrt{n}}\begin{pmatrix}
1 & 0 & \ldots & 0 \\
1 & 1 & \ldots & 0 \\
 \vdots & \vdots & \ddots & \vdots \\
1 & 1 & \ldots & 1 \\
\end{pmatrix}\text{.}$$ PCA corresponds to $A = VD$, where $\Sigma = VD^{2}V^{\top}$ is the singular value decomposition of $\Sigma$. A corresponding decomposition for the BB algorithm is given, for example, in [@br000075].

However, Papageorgiou [@br000100] notes that there are examples where BB and PCA are not giving better results than the forward method in connection with quasi-Monte Carlo. He further shows that the worst case error of integration for a certain class of payoff functions is independent of the path construction. Thus in the sense of worst case error all decompositions are equivalent.

This has been investigated further by Sloan and Wang [@br000120]. The authors show another equivalence principle which roughly states that every decomposition is equally bad and good for QMC, depending on the function that one wants to integrate. For every decomposition $A$ that is good for one payoff function $f$, and every decomposition $\widetilde{A}$ there is another payoff function $\widetilde{f}$ for which $\widetilde{A}$ is equally good.

It is therefore prudent to tailor the decomposition to the problem at hand. This is done, for example, by Imai and Tan [@br000060].

While the possible decompositions of $\Sigma$ provide a clean framework for the study of algorithms for generation of Brownian paths, they are of limited practical value because the matrix-vector multiplication is comparatively slow for all but very small values of $n$. This is the case since general matrix-vector multiplication uses $O\left( n^{2} \right)$ floating point operation (*flops*), while the forward method and the Brownian bridge use only $O\left( n \right)$ flops.

Until recently this has been considered a serious disadvantage of PCA as well, cf. [@br000050]. Yet it has been shown by Scheicher [@br000115], using results from Åkesson and Lehoczky [@br000005], that PCA can be computed using the fast sine transform, thereby using $O\left( n\log\left( n \right) \right)$ flops.

While the importance of the proper choice of the decomposition $AA^{\top} = \Sigma$ for problems arising from quasi-Monte Carlo pricing of financial derivatives is stressed by a number of authors, see e.g. [@br000120; @br000060], there has been a lack of alternatives to the three aforementioned constructions that allow for *fast matrix-vector multiplication*, by which we mean generation of a path using at most $O\left( n\log\left( n \right) \right)$ flops.

The present paper aims to narrow this gap by providing a number of fast constructions of Brownian paths that can be used as alternatives to the three constructions presented above. We know from Sloan and Wang [@br000120] that every one of those constructions will have some payoff functions for which they are especially well suited (and some other payoff functions for which they are ill suited).

The practitioner who is willing to use alternative constructions of Brownian paths therefore is presented with the following alternatives: she might want to tailor a (slow) construction to her special problem in the spirit of Imai and Tan [@br000060], or alternatively, she might want to find out which of the fast constructions presented here is relatively well suited to her problem. Which choice is the better will depend on the special problem. But *if* one of the fast constructions is reasonably close to the optimal construction, its use will reduce the computing time considerably, with essentially the same error.

Papageorgiou [@br000100] observed that there is a one-to-one correspondence between constructions of Brownian paths and orthogonal transforms. We present his simple theorem here since it is essential for all of our constructions.

Theorem 1.1Papageorgiou*If* $\Sigma = SS^{\top}$ *is the Cholesky decomposition of* $\Sigma$*, then any orthogonal transformation* $T$ *on* $\mathbb{R}^{n}$ *defines a decomposition* $\Sigma = ST\left( ST \right)^{\top}$*.Conversely, for every* $n \times n$*-matrix* $A$ *with* $\Sigma = AA^{\top}$*there exists an orthogonal transform* $T$ *such that* $A = ST$*.*

ProofFor any orthogonal transform $T$ we have $T^{\top} = T^{- 1}$, such that $ST\left( ST \right)^{\top} = STT^{\top}S^{\top} = SS^{\top}$.On the other hand, $S$ is invertible, so that for $T = S^{- 1}A$ we have $A = ST$ and $$TT^{\top} = S^{- 1}AA^{\top}\left( S^{- 1} \right)^{\top} = S^{- 1}\Sigma\left( S^{- 1} \right)^{\top} = S^{- 1}SS^{\top}\left( S^{\top} \right)^{- 1} = {id}\text{,}$$ such that $T$ is orthogonal. □

Note that by "orthogonal" we mean that $T$ preserves lengths and angles. Strictly speaking such transforms are "orthonormal", but the term "orthogonal" is more common.

Matrix-vector multiplication with $S$ and $S^{- 1}$ can be done in linear time. Searching for decompositions of $\Sigma$ which admit fast matrix-vector multiplication is therefore equivalent to searching for orthogonal matrices that allow for matrix-vector multiplication using at most $O\left( n\log\left( n \right) \right)$ flops.

In this paper we consider a variety of mostly well-known orthogonal transforms that can be computed using at most $O\left( n\log n \right)$ flops.

Any of those constructions can be modified by combining them with permutations, which do not require any flops at all. Also the multiplication of two or more orthogonal matrices is orthogonal, and the product matrix admits fast matrix-vector multiplications if all of the factors admit this, provided the number of factors is small compared to the number of dimensions.

The remainder of the paper is structured as follows: Section [2](#s000010){ref-type="sec"} reviews the principal component analysis (PCA) construction and the general Brownian bridge (BB) construction, whereby we give some generalizations like interpolation and $m$-step construction, that allow the combination of BB with PCA and other fast constructions.

In the main Section [3](#s000040){ref-type="sec"} we present a number of fast orthogonal transforms derived from fast Fourier transform (FFT) type constructions. A special emphasis is on the fast cosine transform which, as it turns out, corresponds to a construction which is close to PCA. The proof of proximity is one of the main results of our paper.

Section [4](#s000070){ref-type="sec"} contains some fast constructions that do not fit the previous schemes, and connections to previous methods are highlighted.

Some numerical examples are computed in Section [5](#s000090){ref-type="sec"} to illustrate the advantages of some of the proposed methods. Section [6](#s000115){ref-type="sec"} concludes.

Since a large portion of the paper consists of listing fast construction methods, it is probably of interest to the reader to know which of those constructions have previously been known and which are new. This shall be done here, though we need to tread carefully: due to the vast applicability of Brownian motion and simulation we concede that it is too easy to overlook even significant contributions if they are related to a field foreign to the author and maybe even use different terminology.

While the first construction of PCA in Section [2.1](#s000015){ref-type="sec"} is due to Scheicher [@br000115], the second construction using the corresponding orthogonal transform is new. The idea of using the Brownian bridge for interpolation to generalize an equal-step generation method to non-equal steps seems to have been previously unknown, cf. Section [2.3](#s000025){ref-type="sec"}, and the same holds for the general $m$-step forward and BB methods.

To the authors best knowledge, none of the orthogonal transforms presented in Section [3](#s000040){ref-type="sec"} have previously been used for the (fast) construction of Brownian paths, with the sole exception of Scheicher's use of the DST-I to calculate the PCA. But it is worth noting that if one uses DST-I as the underlying orthogonal transform, cf. [Theorem 1.1](#e000005){ref-type="statement"}, then one obtains a construction different from Scheicher's. Nevertheless all of those transforms are of course well known, as are their respective properties. The fact that the construction obtained from DCT-IV is close to PCA is entirely new and consequently so is the idea to use this construction as a faster substitute for PCA when constructing Lévy paths.

The facts from Section [4](#s000070){ref-type="sec"} seem to be folklore to a large extent. However, the construction obtained from the Haar transform is new and it seems that it has not been observed previously that the Kronecker product can be used to combine different (fast) construction methods for Brownian paths.

2. PCA and Brownian bridge {#s000010}
==========================

2.1. PCA construction {#s000015}
---------------------

Åkesson and Lehoczky [@br000005] showed that for $k = 1,\ldots,n$ the $k$-th eigenvalue and eigenvector of the matrix $\Sigma$ are given by $$\lambda_{k} = \left( 4n\sin^{2}\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right) \right)^{- 1}$$ and $v_{k} = \left( v_{k,1},\ldots,v_{k,n} \right)^{\top}$ where $$v_{k,j} = \frac{2}{\sqrt{2n + 1}}\sin\left( \frac{\left( 2k - 1 \right)j}{2n + 1}\pi \right)\text{,}\quad j = 1,\ldots,n\text{,}$$ respectively. Therefore $\Sigma = VD\left( VD \right)^{\top}$, where $V$ is the matrix $$V = \left( v_{1},\ldots,v_{n} \right)\text{,}$$ that is, $V_{jk} = v_{k,j}$, $$\left( Vx \right)_{j} = \sum\limits_{k = 1}^{n}\frac{2}{\sqrt{2n + 1}}\sin\left( \frac{\left( 2k - 1 \right)j}{2n + 1}\pi \right)x_{k}$$ and $D$ is the diagonal matrix that has $\lambda_{1}^{1/2},\ldots,\lambda_{n}^{1/2}$ as its diagonal elements.

It has been observed by Scheicher [@br000115] that $v_{i}$ is essentially the $\left( 2i - 1 \right)$-th basis function of the discrete sine transform (DST-I) in dimension $2n$: Recall the definition of DST-I (provided Section [3.3](#s000055){ref-type="sec"}) on the set $\left\{ 1,\ldots,N \right\}$, $$M_{\text{DST-I}}\left( y \right) ≔ \left( \sqrt{\frac{2}{N + 1}}\sum\limits_{k = 1}^{N}y_{k}\sin\left( \frac{kj}{N + 1}\pi \right) \right)_{j = 1}^{N}\text{.}$$

Therefore $V = \sqrt{2}PM_{\text{DST-I}}Q$, where $P$ is the projection of a $2n$-vector onto its first $n$ elements, $M_{\text{DST-I}}$ is DST-I in dimension $N = 2n$ and $Q$ is the mapping $$Qx ≔ \left( x_{1},0,x_{2},0,\ldots,x_{n},0 \right)\text{,}$$$$\left( M_{\text{DST-I}}Qx \right)_{j} = \sqrt{\frac{2}{N + 1}}\sum\limits_{k = 1}^{n}x_{k}\sin\left( \frac{\left( 2k - 1 \right)j}{N + 1}\pi \right)\text{.}$$ Multiplication with $D,Q$ and $P$ can be done in linear time while multiplication with $M_{\text{DST-I}}$ takes $O\left( 2n\log\left( 2n \right) \right)$ flops (see [@br000130]), therefore the generation of a Brownian path using PCA can be done using $O\left( n\log\left( n \right) \right)$ flops.

Another fast way to construct the PCA is to calculate the corresponding orthogonal transform $T_{PCA} ≔ S^{- 1}DV$ and compute the PCA via $ST_{PCA}$. Theorem 2.1$$T_{PCA} = \sqrt{2}\widehat{P}{\widehat{M}}_{{DCT}\text{-}{III}}\widehat{Q}\text{,}$$*where* ${\widehat{M}}_{{DCT}\text{-}{III}}$ *is the discrete cosine transform* ${DCT}\text{-}{III}$ *in dimension* $2n + 1$*,*$$\widehat{Q}x ≔ \left( 0,x_{1},0,x_{2},0,\ldots,x_{n},0 \right)\text{,}$$*and* $\widehat{P}$ *is the projection of a* $2n + 1$*-vector onto its first* $n$ *elements.*

ProofFor $x \in \mathbb{R}^{n}$$$\left( DVx \right)_{j} = \sum\limits_{k = 1}^{n}\frac{1}{\sqrt{n\left( 2n + 1 \right)}}\left( \sin\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right) \right)^{- 1}\sin\left( \frac{2k - 1}{2n + 1}j\pi \right)x_{k}$$ and for $y \in \mathbb{R}^{n}$$$\left( S^{- 1}y \right)_{j} = \left\{ \begin{matrix}
{\sqrt{n}\, y_{1}} & {j = 1} \\
{\sqrt{n}\left( y_{j} - y_{j - 1} \right)} & {2 \leq j \leq n\text{.}} \\
\end{matrix} \right.$$ Therefore $$\left( S^{- 1}DVx \right)_{j} = \sum\limits_{k = 1}^{n}\frac{1}{\sqrt{2n + 1}}\left( \sin\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right) \right)^{- 1}\left( \sin\left( \frac{2k - 1}{2n + 1}j\pi \right) - \sin\left( \frac{2k - 1}{2n + 1}\left( j - 1 \right)\pi \right) \right)x_{k}\text{.}$$ But using the addition theorem for the sine we get $$\sin\left( \frac{2k - 1}{2n + 1}j\pi \right) - \sin\left( \frac{2k - 1}{2n + 1}\left( j - 1 \right)\pi \right) = \sin\left( \frac{2k - 1}{2n + 1}\left( j - \frac{1}{2} + \frac{1}{2} \right)\pi \right) - \sin\left( \frac{2k - 1}{2n + 1}\left( j - \frac{1}{2} - \frac{1}{2} \right)\pi \right) = 2\cos\left( \frac{2k - 1}{2n + 1}\left( j - \frac{1}{2} \right)\pi \right)\sin\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right)\text{,}$$ such that $$\left( S^{- 1}DVx \right)_{j} = \sqrt{2}\sum\limits_{k = 1}^{n}\sqrt{\frac{2}{2n + 1}}\cos\left( \frac{2k - 1}{2n + 1}\left( j - \frac{1}{2} \right)\pi \right)x_{k}\text{.}$$ Comparing this with the definition in [@br000130] of the DCT-III transform finishes the proof. □

2.2. General Brownian bridge {#s000020}
----------------------------

Suppose $Z = \left( Z_{1},\ldots,Z_{n} \right)$ is a vector of independent standard normal variables.

Consider real numbers $0 = t_{0} < t_{1} < \cdots < t_{n} = 1$. We want to construct a discrete Brownian path $\left( B_{t_{1}},\ldots,B_{t_{n}} \right)$, i.e. a Gaussian vector with covariance matrix $\left( \min\left( t_{j},t_{k} \right) \right)_{j,k}$.

The so-called *forward construction* solves this problem in the following way: $B_{t_{1}} ≔ \sqrt{t_{1}}Z_{1}$ and $B_{t_{k + 1}} ≔ B_{t_{k}} + \sqrt{t_{k + 1} - t_{k}}Z_{k + 1}$ for $k = 1,\ldots,n - 1$. It is easy to check that the vector $\left( B_{t_{1}},\ldots,B_{t_{n}} \right)$ constructed in this way has indeed the required covariance matrix. The construction requires $O\left( n \right)$ flops (note that in a (quasi-)Monte Carlo simulation with many scenarios the values $\sqrt{t_{k + 1} - t_{k}}$, $k = 1,\ldots,n - 1$ need to be computed only once).

An alternative construction is the well-known Brownian bridge, which we will repeat for the convenience of the reader.

Suppose the elements of $\left( B_{t_{1}},\ldots,B_{t_{n}} \right)$ should be computed in the order $B_{t_{\pi{(1)}}},B_{t_{\pi{(2)}}},\ldots,B_{t_{\pi{(n)}}}$ for some permutation $\pi$ of $n$ elements. Consequently, in computing $B_{t_{\pi{(j)}}}$ we need to take into account the previously computed elements. Fortunately at most two of those are of relevance, the one next to $\pi\left( j \right)$ on the left and the one next to $\pi\left( j \right)$ on the right.

Formally define for every $j \in \left\{ 1,\ldots,n \right\}$ two sets, $$L\left( j \right) ≔ \left\{ k:k < \pi\left( j \right)\text{~and~}\pi^{- 1}\left( k \right) < j \right\}$$$$R\left( j \right) ≔ \left\{ k:k > \pi\left( j \right)\text{~and~}\pi^{- 1}\left( k \right) < j \right\}\text{.}$$

That is, $L$ contains all the indices $k$ that are smaller than $\pi\left( j \right)$ and for which $B_{t_{k}}$ has already been constructed and $R$ contains all the indices $k$ that are greater than $\pi\left( j \right)$ and for which $B_{t_{k}}$ has already been constructed. Now define $$l\left( j \right) ≔ \left\{ \begin{matrix}
0 & {\text{if} L_{j} = \operatorname{0\not{}}} \\
{\max L_{j}} & {\text{if} L_{j} \neq \operatorname{0\not{}}} \\
\end{matrix} \right.$$$$r\left( j \right) ≔ \left\{ \begin{matrix}
\infty & {\text{if} R_{j} = \operatorname{0\not{}}} \\
{\min R_{j}} & {\text{if} R_{j} \neq \operatorname{0\not{}}} \\
\end{matrix} \right.$$ and set $B_{t_{0}} = 0$, $$B_{t_{\pi{(j)}}} ≔ \left\{ \begin{matrix}
{B_{t_{l{(j)}}} + \sqrt{t_{\pi{(j)}} - t_{l{(j)}}}Z_{j}} & {\text{if} r\left( j \right) = \infty} \\
{\frac{t_{r{(j)}} - t_{\pi{(j)}}}{t_{r{(j)}} - t_{l{(j)}}}B_{t_{l{(j)}}} + \frac{t_{\pi{(j)}} - t_{l{(j)}}}{t_{r{(j)}} - t_{l{(j)}}}B_{t_{r{(j)}}} + \sqrt{\frac{\left( t_{\pi{(j)}} - t_{l{(j)}} \right)\left( t_{r{(j)}} - t_{\pi{(j)}} \right)}{t_{r{(j)}} - t_{l{(j)}}}}Z_{j}} & {\text{if} r\left( j \right) < \infty\text{.}} \\
\end{matrix} \right.$$ It is straightforward to check that the vector $\left( B_{t_{1}},\ldots,B_{t_{n}} \right)$ constructed in that way has again covariance matrix $\left( \min\left( t_{j},t_{k} \right) \right)_{j,k}$. The functions $l$ and $r$, as well as the factors of $B_{t_{l{(j)}}}$, $B_{t_{r{(j)}}}$, $Z_{j}$, do not depend on the random vector $Z$ so their computation needs to be done only once. In some special cases the functions $l$ and $r$ can be computed explicitly, for example if the $\pi\left( t_{j} \right)$ are the first $n$ elements of the van der Corput sequence or of the $\left\{ k\alpha \right\}$-sequence with $\alpha = \frac{1 + \sqrt{5}}{2}$, see [@br000075]. Since in each step only two of the already constructed values are used it follows that the Brownian bridge construction uses $O\left( n \right)$ floating point operations.

Moreover we see that the forward construction is a special case of the Brownian bridge construction where $\pi\left( j \right) = j$ for all $j$.

The classical Brownian bridge construction as presented in Caflisch and Morokoff [@br000090] corresponds to the setup $n = 2^{L}$, $$t_{k} = k2^{- L}\text{,}\quad k = 0,\ldots,n$$ where the $B_{t_{k}}$ are constructed in the order $B_{1}$, $B_{1/2}$, $B_{1/4}$, $B_{3/4}$, $B_{1/8}$, $B_{3/8}$, $B_{5/8}$, $B_{7/8}$, $\ldots$. This is therefore another example where the functions $l$ and $r$ can easily be computed.

2.3. Interpolation using Brownian bridge {#s000025}
----------------------------------------

The preceding section inspires another class of constructions: we may start with constructing a possibly rough approximation of the discrete path and fill in the gaps using the Brownian bridge.

More concretely, suppose we have nodes $0 < t_{1} < \cdots < t_{N} = 1$ and we want to generate a sample of the corresponding discrete Brownian path, $B_{t_{1}},\ldots,B_{t_{N}}$. We proceed as follows: choose some convenient natural number $n$ of roughly the same or a smaller magnitude as $N$ and use some fast construction to generate $B_{1/n},\ldots,B_{1}$. Then partition the set $\left\{ t_{1},\ldots,t_{N} \right\}$ into subsets $$\left\{ t_{1},\ldots,t_{N} \right\} = T \cup \bigcup\limits_{k = 1}^{n}T_{k}\text{,}$$ where $T_{k} ≔ \left\{ t_{j}:\frac{k - 1}{n} < t_{j} < \frac{k}{n} \right\}$ and $T ≔ \left\{ t_{j}:t_{j} = \frac{k}{n}\;\text{~for~some~}0 < k \leq n \right\}$. The values $B_{t}$ for $t \in T$ are already known. For every $0 < k \leq n$ we now compute $\left\{ B_{t}:t \in T_{k} \right\}$ using the Brownian bridge construction.

The fast construction we used above to generate $B_{1/n},\ldots,B_{1}$ may be any of the constructions presented in Section [3](#s000040){ref-type="sec"} or PCA. Yet it may as well be something entirely different: suppose we have some $n \times n$ matrix $A$ with $AA^{\top} = \Sigma^{(n)}$ and with $n$ being of the order of magnitude of $\sqrt{N}$. Then $\left( B_{1/n},\ldots,B_{1} \right)$ can be constructed in $n^{2} \approx N$ steps (multiplication of an $n$-vector with $A$) and interpolation using the Brownian bridge can be done in at most $O\left( N \right)$ steps.

Another possible solution to the problem of generating a Brownian path with unequally spaced time nodes is presented by Keiner and Waterhouse [@br000070]. They describe an approximate PCA that constructs paths with unequal time-steps.

2.4. $m$-step forward method {#s000030}
----------------------------

Suppose we want to generate a discrete Brownian path on $\left\{ \frac{1}{n},\frac{2}{n},\ldots,1 \right\}$ where $n = n_{1}m$ for some $n_{1},m \in \mathbb{N}$. There is a straightforward generalization of the forward method whereby $m$ forward steps are generated in one "forward leap" in the following way: fix any decomposition $\Sigma^{(m)} = AA^{\top}$. $$\left( B_{{({(k + 1)}n_{1} + 1)}/n},\ldots,B_{{({(k + 1)}n_{1} + m)}/n} \right) ≔ \left( B_{{(kn_{1})}/n},\ldots,B_{{(kn_{1})}/n} \right) + \sqrt{1/n_{1}}AZ_{k + 1}\text{,}$$ where $Z_{1},Z_{2},\ldots$ are independent standard Gaussian vectors of dimension $m$. Clearly the number of flops needed is $n_{1}$ times the number needed by one multiplication by $A$. In any case it is less then $n_{1}m^{2}$ which is $O\left( n \right)$ for constant $m$. If $A$ allows for fast multiplication, the number of operations is $O\left( n\log\left( m \right) \right)$.

2.5. $m$-step Brownian bridge {#s000035}
-----------------------------

Similarly to the $m$-step forward method, we can generalize the Brownian bridge construction to an $m$-step method: suppose we want to generate a discrete Brownian path on $\left\{ \frac{1}{n},\frac{2}{n},\ldots,1 \right\}$ where $n = \left( m + 1 \right)^{\nu}$ for some $m,\nu \in \mathbb{N}$.

Let $\left( B_{\frac{1}{m + 1}}^{0},B_{\frac{2}{m + 1}}^{0},\ldots,B_{1}^{0} \right)$ be a discrete Brownian path and consider the Gaussian vector $$X = \left( X_{\frac{1}{m + 1}},\ldots,X_{\frac{m}{m + 1}} \right) ≔ \left( B_{\frac{1}{m + 1}}^{0} - \frac{1}{m + 1}B_{1}^{0},\ldots,B_{\frac{m}{m + 1}}^{0} - \frac{m}{m + 1}B_{1}^{0} \right)\text{.}$$ It is easy to see that for any standard normal variable $Z$ independent of $B$ the vector $$\left( X_{\frac{1}{m + 1}} + \frac{1}{m + 1}Z,X_{\frac{2}{m + 1}} + \frac{2}{m + 1}Z,\ldots,X_{\frac{m}{m + 1}} + \frac{m}{m + 1}Z,Z \right)$$ is a discrete Brownian path. Moreover the covariance matrix of $\left( X_{\frac{1}{m + 1}},X_{\frac{2}{m + 1}},\ldots,X_{\frac{m}{m + 1}} \right)$ is given by $$\Gamma^{(m)} = \left( \frac{\min\left( j,k \right)}{m + 1} - \frac{jk}{\left( m + 1 \right)^{2}} \right)_{j,k = 1}^{m}\text{.}$$ We call any Gaussian vector $X = \left( X_{\frac{1}{m + 1}},\ldots,X_{\frac{m}{m + 1}} \right)$ with covariance matrix $\Gamma^{(m)}$ a *discrete Brownian bridge* on $\left\{ \frac{1}{m + 1},\ldots,\frac{m}{m + 1} \right\}$. We write $\Gamma$ instead of $\Gamma^{(m)}$ if there is no danger of confusion. $\Gamma$ is a symmetric positive definite matrix which can therefore be decomposed as $\Gamma = CC^{\top}$.

This provides us with a generalization of the classical Brownian bridge algorithm where in every step (except the first one) $m$ points are generated at once.

Let $z_{0}$ be a standard normal variable and let $$\left\{ Z_{j,k}:j = 0,\ldots,\nu - 1,k = 0,\ldots,m^{j} - 1 \right\}$$ be a collection of independent standard normal vectors of dimension $m$.

Define $B_{1} ≔ z_{0}$, subsequently $$\left( B_{m^{- 1}},\ldots,B_{{(m - 1)}m^{- 1}} \right)^{\top} = CZ_{0,0}\text{,}$$ further for $j = 1,\ldots,\nu - 1$$$\left( B_{m^{- j - 1} + km^{- j}},\ldots,B_{{(m - 1)}m^{- j - 1} + km^{- j}} \right)^{\top} = B_{km^{- j}}\left( 1,\ldots,1 \right)^{\top} + \frac{1}{m}\left( B_{{(k + 1)}m^{- j}} - B_{km^{- j}} \right)\left( 1,\ldots,m - 1 \right)^{\top} + \sqrt{m^{- j}}CZ_{j,k}\text{,}$$ for $k \in \left\{ 0,\ldots,m^{j} - 1 \right\}$. It can be shown that this construction generates a Brownian path as required.

Let us now turn to the number of multiplications needed. Multiplication of a vector by $C$ can be done using at most $m^{2}$ scalar multiplications, we need $1 + m + \cdots + m^{\nu - 1}$ such matrix-vector multiplications, giving a total of $m^{2}\frac{m^{\nu} - 1}{m - 1} \approx mn$ scalar multiplications. If matrix-vector multiplication can be done in $m\log\left( m \right)$ steps, then we get $\log\left( m \right)n$. Both results are $O\left( n \right)$ for constant $m$.

The last sentence begs the question under which circumstances matrix-vector multiplication can be done in $m\log\left( m \right)$ steps. We do not know of a direct decomposition of $\Gamma^{(m)}$ that makes that possible, but every fast decomposition of $\Sigma^{(m + 1)}$ gives rise to a fast generation of a discrete Brownian bridge via Eq. [(2)](#fd000150){ref-type="disp-formula"}, using $O\left( \left( m + 1 \right)\log\left( m + 1 \right) \right) = O\left( m\log\left( m \right) \right)$ operations.

A possible drawback is that we use $m + 1$ random numbers to generate an $m$-dimensional Gaussian vector, so that one random number is "lost" or, in other words, the dimension of the integration problem becomes higher with each step.

But this can be remedied to a certain degree: note that if $X$ is constructed from $B^{0}$ according to Eq. [(2)](#fd000150){ref-type="disp-formula"}, then $X$ is uncorrelated and therefore independent of $B_{1}^{0}$: $$\mathbb{E}\left( \left( B_{\frac{k}{m + 1}}^{0} - \frac{k}{m + 1}B_{1}^{0} \right)B_{1}^{0} \right) = \frac{k}{m + 1} - \frac{k}{m + 1} = 0\text{.}$$ Thus $B_{1}^{0}$ may be "recycled" as a component of one of the $Z_{j,k}$ in a later step and so only in the very last step there is one random variable wasted. So if the $m$-step Brownian bridge construction is used to generate a discrete Brownian path on $\left\{ \frac{1}{n},\frac{2}{n},\ldots,1 \right)$ with $n = \left( m + 1 \right)^{\nu}$, then using the above method we need $n + 1 = \left( m + 1 \right)^{\nu} + 1$ random numbers, i.e. the integration problem becomes $\left( n + 1 \right)$-dimensional.

Finally we want to remark that the generation of $X$ above can be used for $m$-step interpolation, analog to the generalized Brownian bridge in Section [2.3](#s000025){ref-type="sec"}.

3. Generation by FFT-type transforms {#s000040}
====================================

Probably the most famous example of a *unitary* transform that allows for fast matrix-vector multiplication is the discrete Fourier transform DFT (with the corresponding fast multiplication algorithm FFT), $$\left( Fx \right)_{j} = \frac{1}{\sqrt{n}}\sum\limits_{k = 1}^{n}x_{k}e^{- {(j - 1)}{(k - 1)}\frac{2\pi i}{n}}\text{,}\quad j = 1,\ldots,n\text{,}$$ for a vector $x \in \mathbb{C}^{n}$. (Usually the indices range from 0 to $n - 1$, but for consistency with other constructions in this paper we chose the above definition.)

There are a many variants of the discrete Fourier transform that map real functions to real functions and are therefore orthogonal. We present some of those variants.

3.1. Modified Fourier transform {#s000045}
-------------------------------

Most Fourier variants that map real functions to real functions have some additional special properties which, for example, make them useful for the fast generation of convolutions. From our point of view we are content with any such variant, indeed we want to have many different variants to choose from for a particular application.

We present one modification of the Fourier transform that maps real functions to real functions. It is well known (and easy to check) that the discrete Fourier transform of a vector $\left( x_{1},\ldots,x_{n} \right)$ is real-valued iff $x_{1} \in \mathbb{R}$ and $x_{n + 2 - k} = {\overline{x}}_{k}$ for $2 \leq k \leq \frac{n}{2} + 1$. Consider the linear map $B$, $$\left( Bx \right)_{k} = \left\{ \begin{matrix}
x_{1} & {k = 1} \\
{x_{k} + ix_{{\lfloor n/2\rfloor} + k}} & {2 \leq k < n/2 + 1} \\
x_{n/2 + 1} & {k = n/2 + 1} \\
{x_{n + 2 - k} - ix_{{\lfloor n/2\rfloor} + n + 2 - k}} & {n/2 + 1 < k \leq n\text{.}} \\
\end{matrix} \right.$$ (Note that one may use this definition for odd $n$ as well.) Then for any standard normal vector $X$, $FBX$ is a (real-valued) standard normal vector.

A drawback of this method is that the FFT algorithm uses complex multiplication, which uses 4 real multiplications. However many fast orthogonal transforms that use only real multiplications have been developed, mainly--as is to be expected--by researchers in the area of signal processing.

3.2. Hartley transform and Hilbert transform {#s000050}
--------------------------------------------

The Hartley transform on $\left\lbrack 0,2\pi \right\rbrack$ is a variant of the Fourier transform where the orthonormal basis is given by the functions $\left. x\mapsto\frac{1}{\sqrt{2\pi}}\left( \cos\left( kx \right) + \sin\left( kx \right) \right) \right.$, $k \in \mathbb{Z}$. See [@br000025; @br000030].

It is easy to see that the Hartley transform on $\left\lbrack 0,2\pi \right\rbrack$ is given by $$\mathcal{H}\left( f \right) ≔ \Re\mathcal{F}\left( f \right) - \Im\mathcal{F}\left( f \right)\text{,}$$ where $\mathcal{F}$ denotes the Fourier transform on $\left\lbrack 0,2\pi \right\rbrack$, $$\mathcal{F}\left( f \right)\left( k \right) = \frac{1}{\sqrt{2\pi}}\int_{0}^{2\pi}f\left( x \right)e^{- 2\pi ikx}dx\text{.}$$

In analogy to the Fourier transform there is a discrete version of the Hartley transform together with a fast algorithm, see [@br000135] (or one might just use the FFT and Eq. [(3)](#fd000195){ref-type="disp-formula"}): for a vector $x \in \mathbb{R}^{n}$ define $$\left( Hx \right)_{j} = \frac{1}{\sqrt{n}}\sum\limits_{k = 1}^{n}x_{k}\left( \cos\left( \left( j - 1 \right)\left( k - 1 \right)\frac{2\pi}{n} \right) + \sin\left( \left( j - 1 \right)\left( k - 1 \right)\frac{2\pi}{n} \right) \right) = \Re\left( Fx \right)\left( j \right) - \Im\left( Fx \right)\left( j \right)\text{,}\quad j = 1,\ldots,n\text{.}$$

We identify the linear map $H$ with the corresponding matrix.

Theorem 3.1*The matrix* $H$ *has the following properties:*1.$H$ *is orthogonal;*2.$H$ *is self-inverse,* $H^{- 1} = H$ *;*3.*for any real vector* $x$*the product* $Hx$ *can be computed in* $O\left( n\log\left( n \right) \right)$ *steps.*

A proof can be found in [@br000030].

Another important property of the Fast Discrete Hartley transform is that it can be computed using only real multiplications and additions. Consult Bracewell [@br000030] for more details on the Hartley transform and its discrete version.

Another example of an orthogonal transform is the so-called *Hilbert transform*. We will not go into the details of the Hilbert transform. The definition and basic properties can be found in [@br000030]. There, a result from Pei and Jaw [@br000105] is also stated which says that the discrete Hilbert transform is of the form $$H \circ \pi \circ H\text{,}$$ where $H$ is the discrete Hartley transform and $\pi$ is the permutation $$\pi\left( j \right) = \left\{ \begin{matrix}
0 & {\text{for} j = 1} \\
{n + 1 - j} & {\text{for} j \neq 1\text{.}} \\
\end{matrix} \right.$$ The discrete Hilbert transform is therefore an orthogonal matrix admitting fast multiplication.

3.3. Sine and cosine transforms {#s000055}
-------------------------------

Two of the most important orthogonal transforms that map real input to real output are the sine and cosine transforms. There are four widely used variants of the sine and cosine transform and there is also the so-called $W$-transform which is similarly defined.

For the convenience of the reader we will recall the definitions of the discrete sine and cosine transforms. Wang [@br000130] or Wickerhauser [@br000135] provide details on fast implementations.

The definitions are not entirely uniform throughout the literature. We have chosen definitions which are in a form that already describes an orthogonal transform. Usually, as for example in [@br000130], the indices range from 0 to $n - 1$, but for consistence with the rest of our paper we will let them range from 1 to $n$.

Define $$\kappa_{A}\left( j \right) ≔ \left\{ \begin{matrix}
1 & {\text{if} j \neq A} \\
{1/\sqrt{2}} & {\text{if} j \in A\text{.}} \\
\end{matrix} \right.$$ All matrices below are in $M_{n,n}\left( \mathbb{R} \right)$ for some $n \geq 2$, so the indices $j,k$ range from 1 to $n$. $$\left( M_{\text{DCT-I}} \right)_{j,k} ≔ \sqrt{\frac{2}{n - 1}}\kappa_{\{ 1,n\}}\left( j \right)\kappa_{\{ 1,n\}}\left( k \right)\cos\left( \frac{\pi}{n - 1}\left( j - 1 \right)\left( k - 1 \right) \right)$$$$\left( M_{\text{DCT-II}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ 1\}}\left( j \right)\cos\left( \frac{\pi}{n}\left( j - 1 \right)\left( k - \frac{1}{2} \right) \right)$$$$\left( M_{\text{DCT-III}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ 1\}}\left( k \right)\cos\left( \frac{\pi}{n}\left( j - \frac{1}{2} \right)\left( k - 1 \right) \right)$$$$\left( M_{\text{DCT-IV}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\cos\left( \frac{\pi}{n}\left( j - \frac{1}{2} \right)\left( k - \frac{1}{2} \right) \right)$$$$\left( M_{\text{DST-I}} \right)_{j,k} ≔ \sqrt{\frac{2}{n + 1}}\sin\left( \frac{\pi}{n + 1}jk \right)$$$$\left( M_{\text{DST-II}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ n\}}\left( j \right)\sin\left( \frac{\pi}{n}j\left( k - \frac{1}{2} \right) \right)$$$$\left( M_{\text{DST-III}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ n\}}\left( k \right)\sin\left( \frac{\pi}{n}\left( j - \frac{1}{2} \right)k \right)$$$$\left( M_{\text{DST-IV}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\sin\left( \frac{\pi}{n}\left( j - \frac{1}{2} \right)\left( k - \frac{1}{2} \right) \right)\text{.}$$ Wang introduced another set of orthogonal transforms which was called discrete $W$ transform, or DWT. The abbreviation "DWT" is nowadays more frequently used for *discrete wavelet transform*. $$\left( M_{\text{W-I}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\sin\left( \frac{\pi}{4} + \frac{2\pi}{n}\left( j - 1 \right)\left( k - 1 \right) \right)$$$$\left( M_{\text{W-II}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ n\}}\left( j \right)\sin\left( \frac{\pi}{4} + \frac{2\pi}{n}\left( j - 1 \right)\left( k - \frac{1}{2} \right) \right)$$$$\left( M_{\text{W-III}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\kappa_{\{ n\}}\left( k \right)\sin\left( \frac{\pi}{4} + \frac{2\pi}{n}\left( j - \frac{1}{2} \right)\left( k - 1 \right) \right)$$$$\left( M_{\text{W-IV}} \right)_{j,k} ≔ \sqrt{\frac{2}{n}}\sin\left( \frac{\pi}{4} + \frac{2\pi}{n}\left( j - \frac{1}{2} \right)\left( k - \frac{1}{2} \right) \right)\text{.}$$

One aspect of the cosine transform is especially worth noting: its corresponding Brownian paths agree to a large extent with those of the principal component analysis[1](#fn000005){ref-type="fn"}.

Let $C ≔ M_{\text{DCT-IV}}$, that is $$C_{jk} = \sqrt{\frac{2}{n}}\cos\left( \frac{\pi}{n}\left( k - \frac{1}{2} \right)\left( j - \frac{1}{2} \right) \right)\text{.}$$ We want to compute $SC$, where $S$ is the scaled summation $$S = \frac{1}{\sqrt{n}}\begin{pmatrix}
1 & 0 & \ldots & 0 \\
1 & 1 & \ldots & 0 \\
 \vdots & \vdots & \ddots & \vdots \\
1 & 1 & \ldots & 1 \\
\end{pmatrix}\text{.}$$$$\left( SC \right)_{lk} = \frac{\sqrt{2}}{n}\sum\limits_{j = 1}^{l}\cos\left( \frac{\pi}{n}\left( k - \frac{1}{2} \right)\left( j - \frac{1}{2} \right) \right) = \frac{\sqrt{2}}{n}\sum\limits_{j = 1}^{l}\cos\left( \frac{\pi}{n}\left( k - \frac{1}{2} \right)\left( j - \frac{1}{2} \right) \right)\overset{( \ast )}{=}\frac{1}{\sqrt{n\left( 2n \right)}}\left( \sin\left( \frac{2k - 1}{2n}\frac{\pi}{2} \right) \right)^{- 1}\sin\left( \frac{2k - 1}{2n}l\pi \right)\text{,}$$ where for $\left( \ast \right)$ one writes $\cos\left( x \right) = \frac{e^{ix} + e^{- ix}}{2}$ and computes the resulting geometric sums.

Recall the PCA construction $$\lambda_{k} = \left( 4n\sin^{2}\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right) \right)^{- 1}$$ and $$V_{l,k} = \frac{2}{\sqrt{2n + 1}}\sin\left( \frac{2k - 1}{2n + 1}l\pi \right)\text{,}\quad l = 1,\ldots,n\text{,}$$$$\left( VD \right)_{lk} = \lambda_{k}^{1/2}V_{k,l} = \frac{1}{\sqrt{n\left( 2n + 1 \right)}}\left( \sin\left( \frac{2k - 1}{2n + 1}\frac{\pi}{2} \right) \right)^{- 1}\sin\left( \frac{2k - 1}{2n + 1}l\pi \right)\text{.}$$ We see that for fixed $l,k$ we have $\lim_{n\rightarrow\infty}\left( \left( SC \right)_{lk} - \left( VD \right)_{lk} \right) = 0$. On its own that does not say a lot about the proximity of paths generated with the two different methods, but [Figs. 1 and 2](#f000005 f000010){ref-type="fig"} illustrate that those paths are rather close indeed, so we want to investigate that topic further.

Consider the expected squared Euclidean norm of the difference of two paths generated by two $n \times n$-matrices $P,Q$ from the same set of independent standard normals $\left( X_{1},\ldots,X_{n} \right)$: $$\mathbb{E}\left( \left\| PX - QX \right\|^{2} \right) = \mathbb{E}\left( \sum\limits_{l = 1}^{n}\left( \left( PX \right)_{l} - \left( QX \right)_{l} \right)^{2} \right) = \mathbb{E}\left( \sum\limits_{l = 1}^{n}\left( \sum\limits_{k = 1}^{n}\left( P - Q \right)_{lk}X_{k} \right)^{2} \right) = \sum\limits_{l = 1}^{n}\sum\limits_{k = 1}^{n}\left( P - Q \right)_{lk}^{2} ≕ d_{n}\left( P,Q \right)^{2}\text{.}$$ The last expression is the square of the Euclidean norm of $P - Q$ in $\mathbb{R}^{n^{2}}$ which is known as the *Hilbert--Schmidt norm* of $P - Q$.

It turns out that $SC$ is close to $VD$ in the Hilbert--Schmidt norm: for small $n$ the value of $d_{n}\left( SC,VD \right)$ can be computed numerically. [Fig. 3](#f000015){ref-type="fig"} shows $d_{n}\left( SC,VD \right)^{2}$ for $n$ ranging from 0 to 50. The graph suggests that $d_{n}\left( SC,VD \right)^{2}$ converges for $\left. n\rightarrow\infty \right.$ and that the limit is close to $\frac{1}{3}$.

In fact we have $d_{n}\left( SC,VD \right) < 1$ for all $n \in \mathbb{N}$ and $$\operatorname{lim\ sup}\limits_{n\rightarrow\infty}d_{n}\left( SC,VD \right)^{2} \leq \frac{2\left( 48 - \pi^{2} \right)}{\left( \pi^{2} - 24 \right)^{2}} = 0.381\ldots\text{.}$$

This will be proved below. As a corollary we note that the average of the variances $\left( \mathbb{E}\left( \left( SCX - VDX \right)_{k}^{2} \right) \right)_{k = 1}^{n}$ tends to 0 as $n$ tends to infinity.

Define $$\delta\left( n,k \right) ≔ \sum\limits_{l = 1}^{n}\left( SC - VD \right)_{lk}^{2}$$ and $$u\left( n,x \right) ≔ \frac{1}{nx^{2}}\left( \frac{9\left( \frac{1}{n} + 1 \right)}{\left( \frac{1}{2n} + 1 \right)^{2}\left( 6 - \left( \frac{1}{2n} + 1 \right)^{2}x^{2} \right)^{2}} + \frac{9}{\left( 6 - x^{2} \right)^{2}} - \frac{\sqrt{n}\sqrt{4n + 2}}{4n + 1} \right)\text{.}$$

Proposition 3.2$u$ *is increasing in the second variable on* $\left\lbrack 0,\frac{\pi}{2} \right\rbrack$*.*

ProofCompute the partial derivative of $u$ w.r.t. $x$ as $$\frac{\partial}{\partial k}u\left( n,k \right) = \frac{1}{nx^{3}}\left( \frac{9\left( \frac{1}{n} + 1 \right)\left( 2n + 1 \right)^{2}x^{2}}{\left( \frac{1}{2n} + 1 \right)^{2}n^{2}\left( 6 - \left( \frac{1}{2n} + 1 \right)^{2}x^{2} \right)^{3}} + \frac{6^{2}x^{2}}{\left( 6 - x^{2} \right)^{3}} - \frac{2 \cdot 9}{\left( 6 - x^{2} \right)^{2}} - \frac{2 \cdot 9\left( \frac{1}{n} + 1 \right)}{\left( \frac{1}{2n} + 1 \right)^{2}\left( 6 - \left( \frac{1}{2n} + 1 \right)^{2}x^{2} \right)^{2}} + \frac{2\sqrt{n}\sqrt{4n + 2}}{4n + 1} \right)\text{.}$$ We need to show that $\frac{\partial}{\partial k}u\left( n,x \right) > 0$, which is equivalent to $$\widehat{u}\left( n,x \right) ≔ nx^{3}\frac{\partial}{\partial k}u\left( n,x \right) > 0\text{.}$$ (Note that $x > 0$). This in turn is shown by proving that $\widehat{u}\left( n,0 \right) > 0$ and $\frac{\partial}{\partial k}\widehat{u}\left( n,0 \right) > 0$. We leave those last calculations to the reader. □

Lemma 3.3*For all* $n$ *and all* $1 \leq k \leq n$ *we have*$$\delta\left( n,k \right) < u\left( n,\frac{2k - 1}{2n + 1}\frac{\pi}{2} \right)\text{.}$$

ProofThe sum $\delta\left( n,k \right) ≔ \sum_{l = 1}^{n}\left( SC - VD \right)_{lk}^{2}$ can be computed by writing $\sin\left( x \right) = \frac{1}{2i}\left( e^{ix} - e^{- ix} \right)$ for $x \in \left\{ \frac{2k - 1}{2n}l\pi,\frac{2k - 1}{2n + 1}l\pi \right\}$ so that the sum becomes a sum of 4 geometric sums. The result can be simplified to $$\delta\left( n,k \right) = \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{8n^{2}} + \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{4n} + \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)}{4n\left( 2n + 1 \right)} + \frac{2\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)}{4\left( 2n + 1 \right)} - \frac{\sqrt{2}\csc\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{4n\left( 2n + 1 \right)} \right)}{4n\sqrt{n\left( 2n + 1 \right)}} - \frac{\sin\left( \frac{\pi\left( 2k - 1 \right)\left( 2n + 1 \right)}{2n} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{2n} \right)\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{8n^{2}} - \frac{\sin\left( \pi\left( 2k - 1 \right) \right)\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{2n + 1} \right)}{4n\left( 2n + 1 \right)} + \frac{\sqrt{2}\sin\left( \frac{\pi\left( 2k - 1 \right)\left( 4n + 1 \right)}{4n} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)\left( 4n + 1 \right)}{4n\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{4n\sqrt{n\left( 2n + 1 \right)}}\text{.}$$ We make use of the facts that $\sin\left( \pi\left( 2k - 1 \right) \right) = 0$, $\sin\left( \frac{\pi\left( 2k - 1 \right)\left( 2n + 1 \right)}{2n} \right) = - \sin\left( \frac{\pi\left( 2k - 1 \right)}{2n} \right)$ and $\sin\left( \frac{\pi\left( 2k - 1 \right)\left( 4n + 1 \right)}{4n} \right) = - \sin\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)$ for integers $k,n$, such that $$\delta\left( n,k \right) = \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{8n^{2}} + \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{4n} + \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)}{4n\left( 2n + 1 \right)} + \frac{2\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)}{4\left( 2n + 1 \right)} - \frac{\sqrt{2}\csc\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)}{4n\left( 2n + 1 \right)} \right)}{4n\sqrt{n\left( 2n + 1 \right)}} + \frac{\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right)}{8n^{2}} - \frac{\sqrt{2}\csc\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)\csc\left( \frac{\pi\left( 2k - 1 \right)\left( 4n + 1 \right)}{4n\left( 2n + 1 \right)} \right)}{4n\sqrt{n\left( 2n + 1 \right)}}\text{,}$$ that is, $$\delta\left( n,k \right) = \frac{1}{4n}\left( \frac{1}{n} + 1 \right)\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{4n} \right) + \frac{1}{4n}\csc^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right) - \frac{\sqrt{2}\csc\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)}{4n\sqrt{n\left( 2n + 1 \right)}}\left( \csc\left( \frac{\pi\left( 2k - 1 \right)}{4n\left( 2n + 1 \right)} \right) + \csc\left( \frac{\pi\left( 2k - 1 \right)\left( 4n + 1 \right)}{4n\left( 2n + 1 \right)} \right) \right)\text{.}$$ Next we use that for $x \in \left( 0,\frac{\pi}{2} \right)$ one has $x - \frac{x^{3}}{6} < \sin\left( x \right) < x$ and therefore $$\frac{1}{x - \frac{x^{3}}{6}} > \csc\left( x \right) > \frac{1}{x}\text{,}$$ such that, after some simplifications we get the estimate $$\delta\left( n,k \right) < \frac{1}{n\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)^{2}}\left( \frac{9\left( \frac{1}{n} + 1 \right)}{\left( \frac{1}{2n} + 1 \right)^{2}\left( 6 - \left( \frac{1}{2n} + 1 \right)^{2}\left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)^{2} \right)^{2}} + \frac{9}{\left( 6 - \left( \frac{\pi\left( 2k - 1 \right)}{2\left( 2n + 1 \right)} \right)^{2} \right)^{2}} - \frac{\sqrt{n}\sqrt{4n + 2}}{4n + 1} \right)\text{.}$$ Therefore $$\delta\left( n,k \right) < u\left( n,\frac{2k - 1}{2n + 1}\frac{\pi}{2} \right)\text{,}$$ □

Theorem 3.4*The sequence* $d_{n}\left( SC,VD \right)^{2}$ *is bounded by* $\frac{2\left( 48 - \pi^{2} \right)}{\left( \pi^{2} - 24 \right)^{2}}$*.*

ProofWe know from [Proposition 3.2](#e000030){ref-type="statement"} that $u$ is increasing in $x$ and from [Lemma 3.3](#e000040){ref-type="statement"} that $\delta\left( n,k \right) < u\left( n,\frac{2k - 1}{2n + 1}\frac{\pi}{2} \right)$ for all $k \in \left\{ 1,\ldots,n \right\}$, so we conclude that $\delta\left( n,k \right) < u\left( n,\frac{\pi}{2} \right)$ and therefore $$d_{n}\left( SC,VD \right)^{2} = \sum\limits_{l = 1}^{n}\sum\limits_{k = 1}^{n}\left( SC - VD \right)_{lk}^{2} = \sum\limits_{k = 1}^{n}\delta\left( n,k \right) < n\, u\left( n,\frac{\pi}{2} \right)\text{.}$$ It is readily checked that $$\lim\limits_{n\rightarrow\infty}n\, u\left( n,\frac{\pi}{2} \right) = \frac{2\left( 48 - \pi^{2} \right)}{\left( \pi^{2} - 24 \right)^{2}} = 0.381\ldots$$ so that $d_{n}\left( SC,VD \right)$ is indeed bounded. □

Some extra work is needed to find a uniform bound on $d_{n}\left( SC,VD \right)$:

Theorem 3.5$d_{n}\left( SC,VD \right) < 1$*.*

ProofObserve that $\frac{\sqrt{4n + 2}\sqrt{n}}{4n + 1} \geq \frac{1}{2} - \frac{1}{64n^{2}}$, so that from [(4)](#fd000350){ref-type="disp-formula"} to [(5)](#fd000400){ref-type="disp-formula"}$$d_{n}\left( SC,VD \right)^{2} < n\, u\left( n,\frac{\pi}{2} \right) \leq \frac{4}{\pi^{2}}\left( \frac{\left( \frac{1}{n} + 1 \right)}{\left( \frac{1}{2n} + 1 \right)^{2}\left( 2 - \left( \frac{1}{2n} + 1 \right)^{2}\frac{\pi^{2}}{12} \right)^{2}} + \frac{1}{\left( 2 - \frac{\pi^{2}}{12} \right)^{2}} - \left( \frac{1}{2} - \frac{1}{64n^{2}} \right) \right) \leq \frac{4}{\pi^{2}}\left( \frac{2}{\left( \frac{1}{2n} + 1 \right)^{2}\left( 2 - \left( \frac{1}{2n} + 1 \right)^{2}\frac{\pi^{2}}{12} \right)^{2}} - \left( \frac{1}{2} - \frac{1}{64n^{2}} \right) \right)\text{.}$$ Now use the estimates $\frac{4}{\pi^{2}} < \frac{41}{100}$ and $\frac{\pi^{2}}{12} < \frac{83}{100}$ to verify that $n\, u\left( n,\frac{\pi}{2} \right) < 1$ for all $n \geq 3$.For $n = 1,2$ direct computation shows that $d_{n}\left( SC,VD \right) < 1$. □

3.4. Intermezzo: simulation of Lévy processes {#s000060}
---------------------------------------------

A *Lévy process* $L$ is a stochastic process with the following properties: 1.$L_{0} = 0$ a.s.;2.for fixed $t > 0$ the random variables $L_{s + t} - L_{s}$ have the same distribution for all $s \geq 0$;3.for $0 \leq t_{1} \leq t_{2} \leq \cdots \leq t_{n}$ the random variables $L_{t_{1}},L_{t_{2}} - L_{t_{1}},\ldots,L_{t_{n}} - L_{t_{n - 1}}$ are independent.4.$L$ is continuous in probability. This is a generalization of Brownian motion, which is obtained if $L_{t}$ has normal distribution with mean 0 and variance $t$.

The distribution of $L_{t}$ for arbitrary $t$ is already determined by the distribution of $L_{1}$: if $\psi$ is the generating function of $L_{1}$, then by definition of a Lévy process and elementary properties of characteristic functions, $L_{m/n}$ has characteristic function $\psi^{m/n}$ for any rational number $m/n > 0$ and by continuity, $L_{t}$ has characteristic function $\psi^{t}$ for any real number $t > 0$. See for example [@br000110].

If we want to know the values of a sample path for $t = \frac{1}{n},\frac{2}{n},\ldots,1$ we can use a direct generalization of the forward method, provided that we know the inverse of the cumulative probability distribution function (CDF) of $L_{1/n}$: let $F$ be the CDF of $L_{1/n}$ and $F^{- 1}$ its inverse. Given independent $\left( 0,1 \right)$-uniform random variables $U_{1},\ldots,U_{n}$ define $$L_{1/n} ≔ F^{- 1}\left( U_{1} \right)$$$$L_{2/n} ≔ F^{- 1}\left( U_{1} \right) + F^{- 1}\left( U_{2} \right)$$$$\vdots$$$$L_{1} ≔ F^{- 1}\left( U_{1} \right) + F^{- 1}\left( U_{2} \right) + \cdots + F^{- 1}\left( U_{n} \right)\text{.}$$ In contrast, the general construction underlying all the generation methods mentioned so far in this paper do not generalize to Lévy processes since those constructions depend heavily on the normality of the generated vector.[2](#fn000010){ref-type="fn"}

It has first been observed in [@br000080] that constructions of Brownian paths like BB and PCA can be generalized to constructions of paths of Lévy processes in the way described below.

Let $\Phi$ denote the standard normal distribution function, i.e. $$\Phi\left( x \right) ≔ \int_{- \infty}^{x}\frac{1}{\sqrt{2\pi}}\exp\left( - \frac{y^{2}}{2} \right)dy\text{.}$$ Given any Brownian path $\left( B_{0},B_{1/n},\ldots,B_{1} \right)$ we have that the random variables ${\widehat{U}}_{k} ≔ \Phi\left( n^{1/2}\left( B_{k/n} - B_{{(k - 1)}/n} \right) \right)$, $k \in \left\{ 1,\ldots,n \right\}$ are independent and uniformly distributed on $\left( 0,1 \right)$. These transformed uniform variables ${\widehat{U}}_{1},\ldots,{\widehat{U}}_{n}$ can now be used in the forward construction of a Lévy path.

Thus every orthogonal transform $T$ can be used for a construction of Lévy paths using the following algorithm: 1.Generate independent, standard normal random variables $X_{1},\ldots,X_{n}$;2.Compute $Y = TX$;3.Compute independent, uniformly distributed random variables ${\widehat{U}}_{1},\ldots,{\widehat{U}}_{n}$ via ${\widehat{U}}_{k} = \Phi\left( Y_{k} \right)$, $k = 1,\ldots,n$;4.Compute a discrete path of $L$ via the forward construction, i.e. $L_{1/n} = F^{- 1}\left( {\widehat{U}}_{1} \right),L_{2/n} = F^{- 1}\left( {\widehat{U}}_{1} \right) + F^{- 1}\left( {\widehat{U}}_{2} \right),\ldots,L_{1} = F^{- 1}\left( {\widehat{U}}_{1} \right) + \cdots + F^{- 1}\left( {\widehat{U}}_{n} \right)$.

If the $T$ in that algorithm corresponds to Brownian path generation using PCA, then the above method can be viewed as a kind of PCA for the Lévy process. In that case $B = VDX$ and ${\widehat{U}}_{k} = \Phi\left( Y_{k} \right)$, where $Y = S^{- 1}VDX$.

If, on the other hand, we use for $T$ the cosine transform of Section [3.3](#s000055){ref-type="sec"}, this will give a result that is close to the one obtained from PCA while using less than half the number of flops, since $PCA$ involves the sine transform on $\mathbb{R}^{2n}$.

This method is actually entirely general: every simulation algorithm can use standard normals instead of uniform random variables, and it is always possible to apply an orthogonal transform to those standard normals. From the probabilistic point of view the resulting integration problem will be equivalent to the original one, but it may still hold that convergence when using QMC or stratification becomes faster.

3.5. Walsh--Hadamard transform {#s000065}
------------------------------

Walsh functions and the corresponding transform were first introduced by Walsh in [@br000125]. Subsequently, Fine [@br000040] augmented the theory and made the connection to classical Fourier theory by observing that the Walsh functions are just the characters of the group $\left( \left\lbrack 0,1 \right), \oplus \right)$, where $\oplus$ is digit-wise addition modulo 2. Like the Fourier transform, the Walsh--Hadamard transform has a discrete version which will be outlined below.

Consider the group structure on $\left( 0,\ldots,2^{n} - 1 \right)$ that is induced by bit-wise operations in base 2. For example $3 \oplus 3 = \left( 1,1 \right)_{2} +_{\mathbb{Z}_{2}^{n}}\left( 1,1 \right)_{2} = \left( 0 \right)_{2} = 0$, $7 \oplus 4 = \left( 1,1,1 \right)_{2} +_{\mathbb{Z}_{2}^{n}}\left( 1,0,0 \right)_{2} = \left( 1,1 \right)_{2} = 3$.

Then the characters of the group are of the form $$\chi_{k}\left( j \right) = \left( - 1 \right)^{k \odot j}\text{,}$$ where $\odot$ denotes bit-wise inner product, i.e. $k \odot j = k_{0}j_{0} + \cdots + k_{m}j_{m}$ for $k = k_{0} + k_{1}2 + \cdots + k_{m}2^{m}$ and $j = j_{0} + j_{1}2 + \cdots + j_{m}2^{m}$. The $\chi_{k}$ are also called *(discrete) Walsh functions*. See [@br000035] for details.

From the character property it follows immediately that $\left( \frac{1}{\sqrt{n}}\chi_{0},\ldots,\frac{1}{\sqrt{n}}\chi_{n - 1} \right)$ form an orthonormal basis of $\mathbb{R}^{n}$ and therefore the coordinate transformation from the canonical basis to Walsh functions is orthogonal. Let $W$ denote the corresponding matrix. Multiplication of a vector with $W$ is called the *discrete Hadamard--Walsh transform*.

Theorem 3.6*The matrix* $W$ *has the following properties:*1.$W$ *is orthogonal;*2.$W$ *is self-inverse,* $W^{- 1} = W$ *;*3.*the entries of* $W$ *are all in* $\left\{ - 1,1 \right\}$ *;*4.*for any real vector* $x$*the product* $Wx$ *can be computed using* $O\left( n\log\left( n \right) \right)$ *flops.*

We omit the easy proof as it is mostly a special case of the Kronecker product of orthogonal matrices which is treated in Section [4.3](#s000085){ref-type="sec"}. A detailed description of the fast algorithm can be found in [@br000045].

There is a nice one-line Mathematica implementation of the fast Walsh transform, taken from [@br000140], which we would like to share with the reader:

A disadvantage of the method is that it relies on $n$ being a power of 2. A big advantage is that most of the occurring floating point multiplications are by 1 or $- 1$, with only $n$ multiplications by $\left( \sqrt{n} \right)^{- 1}$.

4. Further approaches {#s000070}
=====================

4.1. Wavelet analysis {#s000075}
---------------------

Another interesting orthogonal transform is the Haar transform as introduced in [@br000055]. We repeat the basic definitions, loosely following [@br000065]. For a vector $x^{0} ≔ \left( x_{1}^{0},\ldots,x_{n}^{0} \right)$ of length $n = 2^{L}$ define $$x_{k}^{1} ≔ \frac{1}{\sqrt{2}}\left( x_{2k - 1}^{0} + x_{2k}^{0} \right)$$$$d_{k}^{1} ≔ \frac{1}{\sqrt{2}}\left( x_{2k - 1}^{0} - x_{2k}^{0} \right)\text{.}$$ It is not hard to convince oneself that the mapping $\left. x^{0}\mapsto\left( x^{1},d^{1} \right) \right.$ constitutes an orthogonal transform.

Moreover, the construction can be repeated on $x^{1}$ and in general we define $$x_{k}^{j} ≔ \frac{1}{\sqrt{2}}\left( x_{2k - 1}^{j - 1} + x_{2k}^{j - 1} \right)$$$$d_{k}^{j} ≔ \frac{1}{\sqrt{2}}\left( x_{2k - 1}^{j - 1} - x_{2k}^{j - 1} \right)$$ for $k = 1,\ldots,2^{L - j}$, $j = 1,\ldots,L$ and it is again easy to see that the mapping $$\left. x^{0}\mapsto\left( x^{L},d^{L},d^{L - 1},d^{L - 2},\ldots,d^{1} \right) \right.$$ is an orthogonal transform. The inverse can be computed similarly: $$x_{2k - 1}^{j - 1} ≔ \frac{1}{\sqrt{2}}\left( x_{k}^{j} + d_{k}^{j} \right)$$$$x_{2k}^{j - 1} ≔ \frac{1}{\sqrt{2}}\left( x_{k}^{j} - d_{k}^{j} \right)$$ for $k = 1,\ldots,2^{L - j}$, $j = L,\ldots,1$.

Let $H$ denote the Haar transform on $1,\ldots,n$ where $n$ is a power of 2. It is well known (and easy to check) that $SH^{- 1}$, where $S$ is again the normalized summation matrix given in Eq. [(1)](#fd000015){ref-type="disp-formula"}, is exactly the classical Brownian bridge construction. General wavelet analysis may therefore be viewed as a generalization of the Brownian bridge construction.

While the inverse Haar transform does not provide us with a new fast construction method for discrete Brownian paths, the Haar transform $H$ does.

The Haar transform as well as its inverse can be computed using $O\left( n\log\left( n \right) \right)$ flops. See [@br000065] for a proof. If some care is taken the transforms can be done using $O\left( n\log\left( n \right) \right)$ additions/subtractions and $n$ multiplications with powers of $\frac{1}{\sqrt{2}}$.

4.2. Block-diagonal orthogonal matrices {#s000080}
---------------------------------------

The following method can be viewed as a generalization of the $m$-step forward method of Section [2.4](#s000030){ref-type="sec"}: Consider a block-diagonal matrix $$C = \begin{pmatrix}
C_{1} & & & \\
 & C_{2} & & \\
 & & \ddots & \\
 & & & C_{k} \\
\end{pmatrix}\text{,}$$ where $C_{1},\ldots,C_{k}$ are orthogonal matrices and $C_{j}$ has dimension $n_{j} \times n_{j}$, such that $C$ has dimension $n \times n$ with $n = n_{1} + \cdots + n_{k}$.

We have that $C$ admits fast matrix-vector multiplication if every $C_{j}$, $j = 1,\ldots,k$, does, or if $\max_{1 \leq j \leq k}n_{j} \leq K\log\left( n \right)$ for some constant $K$.

4.3. Kronecker product of orthogonal matrices {#s000085}
---------------------------------------------

Recall the Kronecker product[3](#fn000015){ref-type="fn"}$A \otimes B$ of two matrices $A$ and $B$: $$\begin{pmatrix}
A_{11} & \ldots & A_{1m} \\
 \vdots & \ddots & \vdots \\
A_{m1} & \ldots & A_{mm} \\
\end{pmatrix} \otimes \begin{pmatrix}
B_{11} & \ldots & B_{1k} \\
 \vdots & \ddots & \vdots \\
B_{k1} & \ldots & B_{kk} \\
\end{pmatrix} ≔ \begin{pmatrix}
{A_{11}B} & \ldots & {A_{1m}B} \\
 \vdots & \ddots & \vdots \\
{A_{m1}B} & \ldots & {A_{mm}B} \\
\end{pmatrix}\text{.}$$Note that $$\left( A_{1} \otimes A_{2} \right)\left( B_{1} \otimes B_{2} \right) = A_{1}B_{1} \otimes A_{2}B_{2}$$ and $$\left( A_{1} \otimes A_{2} \right)^{\top} = A_{1}^{\top} \otimes A_{2}^{\top}$$ and that the Kronecker product of two identity matrices is an identity matrix. From this it follows that the Kronecker product of orthogonal matrices is an orthogonal matrix.

The significance of the Kronecker product for our considerations is that matrix-vector multiplication with the $mk \times mk$-matrix $A \otimes B$ can be done using less than $\left( mk \right)^{2}$ operations. To that end we use the following algorithm: 1.partition the $mk$-dimensional vector $x$ into $k$ vectors ${\widehat{x}}^{1},\ldots,{\widehat{x}}^{k}$ of dimension $m$, ${\widehat{x}}_{j}^{i} ≔ x_{{(i - 1)}k + j}\,$;2.compute $y_{j} ≔ Bx_{j}$ for $j = 1,\ldots,k$;3.concatenate the $y_{j}$'s to an $mk$-dimensional vector $y$;4.partition $y$ into $m$ vectors ${\widehat{y}}_{1},\ldots,{\widehat{y}}_{m}$ of dimension $k$, ${\widehat{y}}_{j}^{i} ≔ y_{{(j - 1)}m + i}\,$;5.compute $z_{j} ≔ A{\widehat{y}}_{j}$ for $j = 1,\ldots,m$;6.concatenate the $z_{j}$'s to an $mk$-dimensional vector $z$. Then it is easy to check that $z = \left( A \otimes B \right)x$. The algorithm needs $O\left( km^{2} + mk^{2} \right) = O\left( mk\left( m + k \right) \right)$ flops, compared to $O\left( \left( mk \right)^{2} \right)$ flops for classical multiplication of $A \otimes B$ with $x$. If multiplication with $A$ and $B$ is fast, then the complexity further reduces to $O\left( km\left( \log\left( m \right) + \log\left( k \right) \right) \right)$.

Suppose $n = 2^{k}$ and that $$C = A_{1} \otimes A_{2} \otimes \cdots \otimes A_{k}$$ for $k\, 2 \times 2$-matrices $A_{1},\ldots,A_{k}$. Then it follows from induction that the multiplication $Cx$, where $x \in \mathbb{R}^{n}$, can be done in $O\left( 2^{k}k \right) = O\left( n\log\left( n \right) \right)$ operations. An example is the aforementioned Walsh--Hadamard transform which may be written as $W_{2} \otimes \cdots \otimes W_{2}$ ($k$ times), where $$W_{2} = \frac{1}{\sqrt{2}}\left( \begin{array}{cr}
1 & 1 \\
1 & {- 1} \\
\end{array} \right)\text{.}$$This is also an easy proof by induction.

5. Numerical examples {#s000090}
=====================

We consider some numerical examples that illustrate the methods described in the paper.

We consider a European-style option, i.e. one that can be exercised only at the expiry date $T$. We assume that the price process is of the form $$S_{t} = S_{0}\exp\left( \sigma B_{t} + \left( r - \frac{\sigma^{2}}{2} \right)t \right)\text{,}\quad t \in \left\lbrack 0,T \right\rbrack\text{,}$$ where $B$ is a standard Brownian motion, that is, we describe the price process of the share in a Black--Scholes market with interest rate $r$ under the risk-neutral measure.

The arbitrage free price of a derivative $X$ depending only on the path of $S$ up to time $T$ is then given by $$p = e^{- rT}\mathbb{E}\left( X \right)\text{,}$$ provided the expectation is finite. See for example [@br000020] for the general theory of derivative pricing in the Black--Scholes setting.

We limit ourselves to the case where $T = 1$ and where $X$ may only depend on $S_{1/n},S_{2/n},\ldots,S_{1}$ or, equivalently, on $B_{1/n},B_{2/n},\ldots,B_{1}$.

In that case Eq. [(7)](#fd000535){ref-type="disp-formula"} takes on the form $$p = \mathbb{E}\left( f\left( B_{1/n},\ldots,B_{1} \right) \right)\text{,}$$ or, in shortened notation, $p = \mathbb{E}\left( f\left( B \right) \right)$, for some function $f$.

5.1. Asian option {#s000095}
-----------------

A classical example is an average value option, also called Asian option, where $$X = \max\left( \frac{1}{n}\sum\limits_{k = 1}^{n}S_{\frac{k}{n}T} - K,0 \right) = \max\left( \frac{1}{n}\sum\limits_{k = 1}^{n}S_{0}\exp\left( \sigma\sqrt{T}B_{\frac{k}{n}} + \left( r - \frac{\sigma^{2}}{2} \right)\frac{k}{n}T \right) - K,0 \right)\text{.}$$

In [Fig. 4](#f000020){ref-type="fig"}, for example, an Asian option is considered with $n = 64$, $r = 0.045$, $\sigma = 0.3$, $S_{0} = K = 100$, $T = 1$. The option price was evaluated using a quasi-Monte Carlo rule based on a randomly shifted Sobol sequence and the standard normal variates were generated by inversion of the cumulative probability distribution function. Each line shows the $\log_{2}$ of the standard deviation divided by the mean over 64 runs of the rule using a different generation method. Along the $x$-axis we plot the $\log_{2}$ of the number of points used by the QMC-rule.

We see that PCA outperforms the forward construction and DST-III. The inverse Haar transform, which is equivalent to the Brownian bridge, can be seen to be quite close to PCA.

This behavior is rather typical for simple options, e.g. for European style options that depend on the price process of a single share only and in an economically meaningful way, PCA seems to be (almost) a panacea. A remarkable exception is the ratchet option considered by Papageorgiou [@br000100], for which the forward method performs much better.

5.2. Double average option {#s000100}
--------------------------

We give an example of an exotic derivative for which one of the methods proposed in this paper gives a slightly better convergence compared to the standard methods.

A *double average option* is an option for which the payoff is a function of $A - A_{1}$, where $A$ and $A_{1}$ are arithmetic averages of the price process over different intervals. See [@br000145] for details on this kind of Asian option. An example payoff of a double average put is the following: $$X = \max\left( \frac{2}{n}\sum\limits_{k = 1}^{n/2}S_{\frac{k}{n}T} - \frac{2}{n}\sum\limits_{k = n/2 + 1}^{n}S_{\frac{k}{n}T},0 \right)\text{.}$$

[Fig. 5](#f000025){ref-type="fig"} shows a comparison between several generation methods, the three classical methods as well as the one corresponding to the DST-III transform. We see that DST-III consistently outperforms the classical methods, though not by a large margin. But it is worth noting that DST-III is significantly better than the methods using $O\left( n \right)$ flops, Brownian bridge and forward construction.

5.3. Weighted Asian option {#s000105}
--------------------------

It is not hard to find theoretical examples that strongly prefer some non-classical generation method, like it is done by Sloan and Wang [@br000120]. We sketch their idea briefly: suppose the payoff of the option, expressed as a function of the Brownian path, is of the form $f\left( B \right) = g\left( w \cdot B \right)$. Suppose further that $w$ is of the special form $$w^{\top} = \left( 1,0,\ldots,0 \right)T^{- 1}S^{- 1}\text{,}$$ where $T$ is a fixed orthogonal transform and $S$ is scaled summation.

If the Brownian motion is generated via the same orthogonal transform $T$, then $B = STZ$, where $Z$ is a vector of independent standard normals, and so $$w \cdot \left( STZ \right) = w^{\top}STZ = \left( 1,0,\ldots,0 \right)T^{- 1}S^{- 1}STZ = \left( 1,0,\ldots,0 \right)Z = Z_{1}\text{,}$$ such that $f\left( B \right) = g\left( Z_{1} \right)$, and the integration problem becomes one-dimensional.

Therefore one may expect a much faster convergence of the QMC rule than for the $n$-dimensional original problem that is obtained for most other generation methods of the Brownian path.

But the same effect may be observed in the (practically more interesting) case that the payoff is only approximately of the form $f\left( B \right) = g\left( w \cdot B \right)$.

Consider a weighted version of the Asian option, that is $$X = \max\left( \sum\limits_{k = 1}^{n}w_{k}S_{\frac{k}{n}T} - K,0 \right) = \max\left( \sum\limits_{k = 1}^{n}w_{k}S_{0}\exp\left( \sigma\sqrt{T}B_{\frac{k}{n}} + \left( r - \frac{\sigma^{2}}{2} \right)\frac{k}{n}T \right) - K,0 \right)\text{,}$$ where $w^{\top} = \left( 1,0,\ldots,0 \right)T_{OTP}^{- 1}S^{- 1}$, and where $T_{OTP}$ is an orthogonal Kronecker product as described in Section [4.3](#s000085){ref-type="sec"}, $T = A \otimes \cdots \otimes A$ (6 times) with $$A = \left( \begin{array}{rr}
{\cos\left( \phi \right)} & {\sin\left( \phi \right)} \\
{- \sin\left( \phi \right)} & {\cos\left( \phi \right)} \\
\end{array} \right)$$and $\phi = 2\pi/3$.[4](#fn000020){ref-type="fn"}[Fig. 6](#f000030){ref-type="fig"} shows the corresponding graphs. We observe that now OTP is the best method by a large margin, while PCA and BB perform roughly as well as the forward method.

[Fig. 7](#f000035){ref-type="fig"} shows the weight $w$ that is used in the payoff function of the weighted Asian option in this particular example. It shows that the option corresponding to this weight is rather exotic indeed.

Remark 5.1As can already be seen from the construction of $w$, there will in general be many orthogonal transforms that make the integration problem (roughly) one-dimensional.

5.4. Weighted Asian option in the NIG Lévy model {#s000110}
------------------------------------------------

We conclude our study of Asian options with considering the same example of a (weighted) Asian option in an NIG Lévy setup. Here the log-increments of the price process have NIG distribution instead of normal distribution, i.e. model [(6)](#fd000530){ref-type="disp-formula"} is replaced by $$S_{t} = S_{0}\exp\left( L_{t} + \left( r - \xi \right)t \right)\text{,}$$ where $L$ is an NIG Lévy process, that is, $L_{t + \Delta t} - L_{t}$ has NIG density $$f_{NIG}\left( x;\alpha,\beta,\delta,\mu \right) = \frac{\alpha}{\pi}\exp\left( \delta\sqrt{\alpha^{2} - \beta^{2}} + \beta\left( x - \mu \right) \right)\frac{\delta K_{1}\left( \alpha\sqrt{\delta^{2} + \left( x - \mu \right)^{2}} \right)}{\sqrt{\delta^{2} + \left( x - \mu \right)^{2}}}\text{,}$$ where $K_{1}\left( x \right)$ is the modified Bessel function of the second kind with index 1, that is $$K_{1}\left( x \right) = \frac{1}{2}\int_{0}^{\infty}\exp\left( - \frac{1}{2}x\left( z + z^{- 1} \right) \right)dz$$ and $\alpha,\beta,\delta,\mu$ are the parameters of the distribution. The parameters $\alpha$ and $\beta$ need to satisfy $0 \leq \left| \beta \right| \leq \alpha$ and $\delta > 0$. $\xi$ is chosen so that $e^{- rT}\mathbb{E}\left( S_{T} \right) = S_{0}$.

The parameters in our numerical example are $\alpha = 18.3$, $\beta = - 1.06$, $\delta = 0.0184$ and $\mu = 0.000434$. Those have been found from the stock price process of General Electrics using maximum likelihood estimation of the data from January 2000 to September 2010. The popularity of Lévy processes derives from the empirical fact that they give a much better fit to market data than the Gaussian Black--Scholes model. See [Fig. 8](#f000040){ref-type="fig"} for a comparison of the respective maximum likelihood estimated densities with observed market frequencies.

[Fig. 9](#f000045){ref-type="fig"} shows the result for the same weighted Asian option as before, but now the NIG model has been used to generate the price process. We see that the behavior is similar to the Black--Scholes case, but it is much less expressed.

6. Conclusions and open questions {#s000115}
=================================

We have given a number of constructions of discrete Brownian paths that provide alternatives to the classical constructions, that is, forward, Brownian bridge and PCA construction. All the constructions presented have the desirable property that they have computational complexity $O\left( n\log\left( n \right) \right)$ or $O\left( n \right)$.

We have used the orthogonal representation of path construction to derive a new method for efficient computation of the PCA construction and an even faster method for approximate computation of the PCA construction.

We provided numerical examples illustrating two main points: (1) there are cases where alternative constructions are more efficient then the classical ones, (2) however for most practical cases PCA is best with the computationally less complex Brownian bridge being close.

We need to stress that the examples presented all depended on only one Brownian path. Imai and Tan [@br000060] provide examples where a generation method other than PCA is far more effective, however their method uses full matrix multiplication. It would be useful to try and find a fast version of their method or, more generally, find a generic algorithm which combines the orthogonal transforms presented here in an optimal (or near optimal) way for a given integration problem.

The author would like to thank Gerhard Larcher, Fritz Pillichshammer and the anonymous referees for valuable comments.

The author is partially supported by the Austrian Science Foundation (FWF), Project P21196.

This is not a complete surprise, since the scaled summation $S$ serves as a kind of integral and the indeterminate integral of the cosine is the sine. Nevertheless the above assertion is not trivial to prove, and it does not hold for DCT-I and DCT-III.

Nevertheless there are Bridge constructions for Lévy processes for which the conditional distribution of $L_{t/2}$ given $L_{t}$ can be obtained explicitly, like for the variance-gamma process, cf. [@br000015].

We confine ourselves to quadratic matrices.

This kind of transform has been dubbed "CRAFOT"--constant rotation angle fast orthogonal transform--by Misans and Terauds [@br000085].

![Two Brownian paths with $n = 2^{6}$, one generated via PCA and one via DCT-IV.](gr1){#f000005}

![Two Brownian paths with $n = 2^{8}$, one generated via PCA and one via DCT-IV.](gr2){#f000010}

![The squared distance $d_{n}\left( SC,V\, D \right)^{2}$ between PCA and DCT-method for $n = 1,\ldots,50$.](gr3){#f000015}

![Asian option.](gr4){#f000020}

![Double average put option.](gr5){#f000025}

![Weighted Asian option.](gr6){#f000030}

![Weight of the Asian option.](gr7){#f000035}

![Fit of NIG and Gaussian distribution to market log-returns.](gr8){#f000040}

![Weighted Asian option in the NIG Lévy model.](gr9){#f000045}
