# Simulating a process

Let $$W$$ be a random variable valued in $$L^2[0,1]$$ (an infinite dimensional function space). Take $$W=\{W(t), t\in[0,1]\}$$ on $$[0,1]$$.

$$W(t)=\sum_{i=1}^\infty e_i(t) N_i, \quad \forall t\in [0,1]$$ where for all $$i\in\mathbb N$$, $$e_i(t)=\sqrt{2}\sin{(i-1/2)\pi t},$$ $$N_i \sim N(0,\lambda_i), \quad \lambda_i=[(i-1/2)\pi]^{-2}, \{N_i\} \text{ independent}.$$

QUESTION How can I generate each process? Each value $$W(t)$$ is the limit of a series. Should I truncate the indices of the summation to $$\{1,\dotsc, k\}$$ where $$k$$ is a fixed "large" integer?

This process is well known and there should be a simple way to perform simulations.

• There must be typographical errors in your expressions, because $\lambda_i$ never appears. As far as your question goes, are you asking about how to understand the theorem (which explicitly tells you how to approximate the process) or how to generate a realization of $W$ (which usually does not invoke the theorem)?
– whuber
Mar 22, 2023 at 12:57
• @whuber $\lambda_i$ is the variance of the gaussian random variable $N_i$. I want to generate a realization of $W$. Mar 22, 2023 at 18:04
• Thank you for pointing out where $\lambda_i$ appears.
– whuber
Mar 22, 2023 at 18:56

## 2 Answers

One approach is to use the theorem. It is an orthogonal expansion designed to give the best possible approximation (in a least squares sense) whenever you cut it off at a finite endpoint.

Notice that since

$$1 = \operatorname{Var}(W(1)) = \sum_{i=1}^\infty e_i(1)^2\operatorname{Var}(N_i) = \sum_{i=1}^\infty 2 \lambda_i,$$

the $$\lambda_i$$ sum to $$1/2.$$ (You can verify this analytically by relating the sum of the $$\lambda_i$$ to $$\zeta(2) = \pi^2/6.$$)

Thus, to approximate the process with a least squares error smaller than $$\epsilon,$$ find the first $$n$$ for which

$$1 - \epsilon \le 2\sum_{i=1}^n \lambda_i,$$

generate $$n$$ iid Normal values $$Z_1,Z_2,\ldots, Z_n,$$ and set $$N_i = \sqrt{\lambda_i} Z_i.$$ Roughly, to improve the precision from $$\epsilon$$ to $$\epsilon/C$$ (with $$C\ge 1$$) you will need $$C$$ times as many terms. (An excellent approximation for $$n$$ is $$0.2/\epsilon.$$)

Here is an example of one realization of this process using a common sequence of realized values of $$N_i$$ with varying levels of precision:

The advantages of this approach are clear: with limited computing time you can generate an approximation of the walk and, given additional time, you can dynamically improve it to supply detail, as illustrated by the progression in this illustration. Moreover, you can zoom in at will merely by keeping the same realized values of the $$N_i$$ and computing the sine functions within arbitrarily narrow intervals. (The further you zoom in, though, the more accurate your approximation will need to be.)

The code in R is particularly simple once you have computed the $$\lambda_i:$$ one line to generate the $$N_i$$ and another line to compute the sum. See the block of code below following "One realization of the process." It finds $$n$$ by a simple search, doubling the search length until $$n$$ is found or the problem is getting too large. Generally you would be fine just rounding $$0.2/\epsilon$$ up to an integer.

n.max <- 1e4
par(mfrow = c(2,2))
for (eps in c(1e-1, 1e-2, 1e-3, 1e-4, 1e-5)) {
#
# Search for n.
#
n <- 1
lambda <- 8/pi^2
remainder <- 1 - lambda
while(remainder > eps && n < n.max) {
i <- seq_len(n) + n
n <- n + length(i)
l <- 2/(pi * (i - 1/2))^2
lambda <- c(lambda, l)
remainder <- remainder - sum(l)
}
if (n > n.max) {
warning("Intended precision not attained with ", n.max, " terms.")
break
}
Lambda <- cumsum(lambda)
n <- which.max(Lambda >= 1 - eps)
lambda <- lambda[1:n]
Lambda <- Lambda[1:n]
#
# One realization of the process.
#
set.seed(17)
N <- rnorm(n) * sqrt(lambda)
W <- function(t, N) {
sqrt(2) * colSums(outer(seq_along(N), t, \(i,x) sin((i - 1/2) * pi * x) * N))
}
curve(W(t, N), 0, 1, xname = "t", n = 1201,
main = bquote(paste(epsilon == .(eps), ", ", n == .(n))))
}
par(mfrow = c(1,1))

• Impressive answer! Thank you. Mar 22, 2023 at 19:58

While not based on the KL theorem you state, Wiener processes can be approximated/simulated as follows, using that the Wiener process is the continuous-time limit of a (scaled) random walk

n <- 5000
u <- rnorm(n)
W <- 1/sqrt(n)*cumsum(u)