Simulating a stochastic integral

Question

I am trying to solve exercise 3.9.10 on p. 66 of Ubbo F. Wiersema's "Brownian Motion Calculus" (John Wiley & Sons, 2008), which asks to simulate the stochastic integral $$ \int_0^1 B(t)\ dB(t) $$ by initially using a partition of $[0, 1]$ into $n = 2^8$ sub-intervals and running $1000$ simulations of the discrete stochastic integral $$ I^{(n)} = \sum_{i = 0}^{n - 1} B(t_i)\left(B(t_{i + 1}) - B(t_i)\right) $$ for this $n$, then repeating this procedure by repeatedly doubling $n$ to $2^{11}$. The exercise asks to compare the results against the mean and variance of the closed form of the integral, namely $$ \int_0^1 B(t)\ dB(t) = \frac{1}{2} B(1)^2 - \frac{1}{2} $$ which are $0$ and $\frac{1}{2}$, respectively (since $B(1)$ is a standard normal random variable).

I have written a simulation in R, as shown below. However, the variances resulting from this simulation are

> v_sim
        8         9        10        11 
0.4771895 0.4304475 0.5260542 0.4664552

which don't appear to converge to the anticipated value $0.5$. The same phenomenon is observed when changing the seed to $2$ and $3$ as well as when the number of sub-intervals is increased to $2^{12}$ to $2^{15}$. What am I doing wrong?

NSIMS <- 1000L # number of simulations of the integral for every "level"
LEVELS <- 8L:11L # level n corresponds to 2^n sample points between 0 and 1
SEED <- 1

set.seed(SEED)

sims <- data.frame(simid = 1L:NSIMS) # sims contains one colum per every level.
                                        # The column simid is a dummy column
                                        # intended to avoid an error message
                                        # when columns are added to sims
                                        # on the fly inside the following loop.
                                        # This column is deleted after the loop
                                        # ends.
for (n in LEVELS)
{
        nticks <- 2^n # nticks is the number of sample points between 0 and 1
        delta <- 1/nticks
        ticks <- seq(from = 0, to = 1, by = delta)
        std <- sqrt(delta)
        sim <- vector(mode = "numeric", length = NSIMS)
        for (j in 1L:NSIMS)
        {
                b <- cumsum(c(0, rnorm(nticks, sd = std))) # b is a simulated
                                                        # Brownian motion path
                                                        # sampled at the tick
                                                        # marks.
                integral <- 0
                for (i in 1L:(length(b) - 1L))
                {
                        integral <- integral + b[i]*(b[i + 1] - b[i])
                }
                sim[j] <- integral
        }
        sims[, as.character(n)] <- sim
}
sims$simid <- NULL

m_sim <- sapply(sims, mean)
v_sim <- sapply(sims, var)

Answer 1

You have two sources of divergence from 0.5, and you're only making changes to one of them; as a result you won't expect to see the sort of convergence you suggest. [Instead you'd expect to see a slight improvement at the smallest $n$ values which then disappears into noise as you increase $n$ further.]

The first source of divergence is recognized in your question -- the approximation of the continuous process by a discrete approximation.

The second source is the random variation due to simulation. At some finite number of simulations ($b$), you expect that the sample variation calculation will be different from the underlying "true" value you'd get with any given discrete approximation as $b\to\infty$.

To see whether this is indeed the issue, you might try a sequence of $n$ values, say n=100,1000,10000 (and perhaps more if you can) at each $n$. Alternatively, you can compute a theoretical variance of the quantities you print out (and hence standard error), from which you can tell whether your results are clearly inconsistent with 1/2 (as $n$ becomes large you should find the inconsistency becomes harder to see because it would require much larger $b$ to identify).

Answer 2

Let $X_i := B(t_i) [B(t_{i+1}) - B(t_i)]$ for $i=0$ to $n-1$. The rv.s $X_i$ are gaussian and independent due to the independence of the increments of $B(t)$. Moreover $\mathbb{E}(X_i) = 0$ and $$ \text{Var}(X_i) = E[X_i^2] = \mathbb{E}[B(t_i)^2] \, \mathbb{E}[(B(t_{i+1}) - B(t_i))^2] = t_i \, (t_{i+1} - t_i), %% = \frac{i}{n} \frac{1}{n} = \frac{i}{n^2} $$ so $\text{Var}(X_i) = i / n^2$ if the points $t_i := i /n$ are used. Then $$ \text{Var}(I^{(n)}) = \sum_{i=0}^{n-1} \text{Var}(X_i) = \frac{1}{n^2} \sum_{i=0}^{n-1} i = \frac{1}{n^2} \, \frac{(n-1)n}{2} = \left[1 - \frac{1}{n}\right] \frac{1}{2}. $$ Up to the term $1/n$, the distribution of $I^{(n)}$ remains nearly the same for all $n$. By drawing a number, say $N_{\text{sim}}$, of independent $I^{(n)}$, we get a sample of size $N_{\text{sim}}$ from a distribution depending only slightly of $n$ and which is nearly the normal with mean $0$ and variance $1/2$. So in the program, the elements of v_sim corresponding to n are (unless n is small) simply the variances of independent samples of size $N_{\text{sim}}$ from $\text{Norm}(0, \,1/2)$.

The number $n$ of points is nearly irrelevant here because of the specific integrand considered. The situation can be compared to the evaluation of the non-stochastic integral $\int_0^1 t \,\text{d}t$ using the trapezoidal rule. The result will not change when more trapezes are used. Obviously things would be different if a different adapted process or function was used as integrand, a denser design $[t_i]_i$ leading then to a smaller bias.

Answer 3

To simulate the convergence, simulate the mean of the squared DIFFERENCE between the integral based on n steps and the integral based on 2n steps. Then the same for the difference between the integrals for 2n steps and 4n steps. Then the difference between 4n steps and 8n steps. Use the SAME Brownian motion trajectory (in the unit time period) for all. The same convergence pattern can be seen by using the variance of the DIFFERENCE between the pairs of integrals mentioned above. Section 3.4 of my book covers this. Ubbo

Answer 4

3.9.10 BdB.R

integral B.dB

rephrase exercise 3.9.10 as follows:

delete line 5 and onward until line 8 sentence that starts with Repeat

also delete sentence starting with Compare

set.seed(123) # optional; positive integer sims <- 1000 # number of simulations

BdB.sim <- function(n) { B.steps <- 2*n # number of time steps in Brownian Motion path

# integrals B.dB based on successively halved time partitions
I.n <-  numeric(sims)  # n steps
I.2n <- numeric(sims)  # 2*n steps

dt <- 1/B.steps
sqrt.dt <- sqrt(dt)
k.max <- 1+B.steps     # last timepoint index of B

# timepoints
# ---------------------------------------
k.seq.n <- seq(from=3,to=k.max,by=2)
length(k.seq.n)

k.seq.2n <- seq(from=2,to=k.max,by=1)
length(k.seq.2n)

k.seq.step <- seq(from=3,to=k.max-1,by=2)
length(k.seq.step)
# ---------------------------------------

for(sim.nr in 1:sims) {
    # each simulation starts by generating a Brownian Motion path of n timesteps, denoted B
    B <- c(0,cumsum(rnorm(B.steps,mean=0,sd=sqrt.dt)))  # Brownian Motion path timestep dt; B[1]=0
    for(k in k.seq.n) {I.n[sim.nr] <- I.n[sim.nr] + B[k-2]*(B[k]-B[k-2])}
    for(k in k.seq.2n) {I.2n[sim.nr] <- I.2n[sim.nr]  + B[k-1]*(B[k]-B[k-1])}
}

mean.sq.diff <- mean((I.n-I.2n)^2)
exact <- 1/(4*n)

return(data.frame(n=n,double=2*n,mean.sq.diff=mean.sq.diff,exact=exact))

}

(n.seq <- c(2^8, 2^9, 2^10,2^11))

trail <- data.frame(n=NA,double=NA,mean.sq.diff=NA,exact=NA)

for(m in 1:length(n.seq)) { out <- BdB.sim(n.seq[m]) trail <- rbind(trail,out) }

trail[-1,] # shows convergence

write.csv(trail[-1,],"3.9.10 BdB convergence.csv",row.names=FALSE) # optional

--------------------------

Ubbo Wiersema 26 Feb 2018