Solution (source code)

= Solution

In the <Lévy construction of Brownian motion>, start with $B_0=0$ and $B_1=Z_0$, where $Z_0$ has the <standard normal distribution>. Suppose values are defined at the endpoints $u=k2^{-n}$ and $v=(k+1)2^{-n}$. Set
$$
B_{(u+v)/2}=\frac{B_u+B_v}{2}+\frac{\sqrt{v-u}}2Z_{n,k},
$$
using fresh <independent> <standard normal random variables> $Z_{n,k}$, also <independent> of everything previously used. Join the values by straight segments to obtain a continuous piecewise linear process $B^{(n)}$ on the dyadic grid.

To verify the <Brownian increments>, write $D=B_v-B_u$ and $\ell=v-u$. The two new increments are $D/2+\sqrt\ell Z_{n,k}/2$ and $D/2-\sqrt\ell Z_{n,k}/2$. If $D\sim N(0,\ell)$, their <variances> are $\ell/2$ and their <covariance> is zero. They are <Jointly Gaussian>, so they are <independent>. Induction proves the required <Brownian increments> on every dyadic grid.

For convergence, successive interpolants differ by tents of disjoint interiors, so
$$
\|B^{(n+1)}-B^{(n)}\|_\infty
=2^{-n/2-1}\max_{0\leq k<2^n}|Z_{n,k}|.
$$
The <Gaussian tail bound> and the <union bound> give
$$
\mathbb P\!\left(\max_{k<2^n}|Z_{n,k}|>n\right)
\leq2^{n+1}e^{-n^2/2}.
$$
These <probabilities> are summable. The <First Borel-Cantelli lemma> and $\sum_n n2^{-n/2}<\infty$ therefore prove <uniform convergence> of the interpolants <almost surely>. Their limit has continuous <sample paths>. Approximating any finite set of times by dyadic times passes the <Jointly Gaussian> increment laws to the limit, giving a standard <Brownian motion> on $[0,1]$.

Equivalently, with $\psi_{n,k}(t)=\int_0^tH_{n,k}(s)\,ds$, these tents give the <Lévy-Ciesielski expansion>
$$
\boxed{B_t=Z_0t+\sum_{n=0}^\infty\sum_{k=0}^{2^n-1}Z_{n,k}\psi_{n,k}(t),
\quad\text{uniformly almost surely on }[0,1].}
$$