How far away is $\max_{x: x \in \{0, \ldots, N\}} |W(x/N)|$ from $\max_{0 \leq t \leq 1} |W(t)|$ ($W(t)$ a Wiener process)? How far away is
$$\max_{x: x \in \{0, \ldots, N\}} \left|W\left(\frac{x}{N}\right)\right|$$
from
$$\max_{0 \leq t \leq 1} |W(t)|$$
In other words, if you simulate a Wiener process over a finite set of equidistant rational numbers in $[0,1]$ and take the maximum of the Wiener process at those points, how far away is that maximum away from the true maximum of the Wiener process (in any sense you want, be it probability, distribution, expectation, almost surely, etc.)?
This matters because I need to simulate this maximum and since I cannot do so over a continuous set on a computer I would like to know how much of an error I make by working with the next-best thing: a very fine discretization of $[0,1]$.
Another approach to simulating the maximum would be to find the maximum over the discrete set, take the two closest neighboring points, then simulate Brownian bridges connecting the points and computing a new maximum. This seems like a lot of computational work, though, and also seems to just move the problem elsewhere.
Any results appreciated.
 A: The convergence of the discretized version $\max_{x \in \{0, \ldots, N\}} |W(x/N)|$ of $M:=\max_{0 \le t \le 1} |W(t)|$ to $M$ will be very slow -- at the rate of $1/\sqrt N$, according to Korolyuk 1961 and Nagaev 1970 (Korolyuk 1961 apparently exists only in Russian, but can be rather easily read using e.g. Google Translate), and then you will have to simulate a very large number $N$ of normal random variables -- just to get one realization of $M$. 
A much more efficient way is to use the following explicit expression for the cdf of $M$: 
\begin{equation}
 P(M\le x)=\frac4\pi\, \sum _{k=0}^{\infty } \frac{(-1)^k}{2 k+1}\,
  \exp \left(-\frac{(2 k+1)^2 \pi ^2}{8
   x^2}\right) \tag{1}
\end{equation}
for $x>0$; see e.g. page 3. The series in (1) converges very fast unless $x$ is very large. On the other hand, if $x$ is large, then $P(M\le x)$ is very close to $1$. 
More specifically, if $x>x_0$, where $x_0=5.8471\ldots$ is the root of the equation $4(1-\Phi(x_0))=1/10^8$ and $\Phi$ is the standard normal cdf, then
\begin{equation}
 0<1-P(M\le x)=P(M>x)\le P(M^+>x)+P(M^->x)=2P(M^+>x)=4(1-\Phi(x))<1/10^8, 
\end{equation}
where $M^\pm:=\max_{0\le t\le1}(\pm W(t))$; the last displayed equality follows immediately 
from another formula on the same page 3. 
On the other hand, for $x\in(0,x_0]$, retaining only the first $11$ summands of the alternating series in (1), we will get the approximation 
\begin{equation}
 P(M\le x)\approx F(x):=\frac4\pi\, \sum _{k=0}^{10} \frac{(-1)^k}{2 k+1}\,
  \exp \left(-\frac{(2 k+1)^2 \pi ^2}{8
   x^2}\right) 
\end{equation} 
with an error less in absolute value than 
\begin{equation}
 \frac4\pi\,\frac1{2\times11+1}\,
  \exp \Big(-\frac{(2\times11+1)^2 \pi ^2}{8x_0^2}\Big)<3/10^{10}. 
\end{equation}
Now, using the inverse transform sampling method, we can get very close approximations $m_1,m_2,\dots$ to iid realizations of $M$ as solutions of the equations $F(m_i)=u_i$ for $m_i$, where $u_1,u_2,\dots$ are iid realizations of a random variable uniformy distributed between $0$ and $1$. It then takes under a second in Mathematica to simulate $1000$ very close approximations $m_1,\dots,m_{1000}$ to iid realizations of $M$: 

