Simulating the Squared Norm of Independent Brownian Motions

Summary

The document considers a process formed by summing the squares of several Brownian motions and derives its stochastic differential using Itô's formula. It then presents a simulation that accumulates the Brownian increments and the corresponding drift, asking why the simulated mean is about three times the expected value.

The key modeling point is that the squared norm of an n-dimensional Brownian motion is a squared Bessel process, whose dynamics include a state-dependent diffusion term. In the posted code, the increment for each squared Brownian component is approximated by twice the current Brownian value times its increment, plus a drift contribution. However, the implementation's Brownian path and increments are misaligned: the path is cumulatively summed from an array whose first column is forced to zero, and the returned path and increment slices do not pair each value with the increment that advances it. The document poses this discrepancy but provides no answer or numerical validation, so it does not establish the precise source of the reported factor or discuss discretization error.

Key ideas

• The sum of squared independent Brownian motions is a squared Bessel process.
• Itô's formula gives the process a constant drift equal to the number of dimensions.
• The diffusion increment depends on the current Brownian values, so paths and increments must be aligned in simulation.
• The document reports a mean discrepancy but does not resolve it or validate the proposed implementation.

Tags

# Simulating sum of squared brownian motions process


# Simulating sum of squared brownian motions process












I'm trying to simulate the following stochastic process:

\begin{equation} R_t = \sum_{i=1}^nB_{i,t}^2 \end{equation}

which has the following dynamics:

\begin{equation} \begin{aligned} dR_t = \sum_{i=1}^n 2B_{i,t}dB_{i,t} + \frac{1}{2}\sum_{i=1}^n 2 dt \\ = ndt + \sum_{i=1}^n 2B_{i,t}dB_{i,t} \end{aligned} \end{equation}

for this I created a function that simulates the brownian motion process and then calculated the dymacis of the process:

```
def brownian_motion(T: int, steps: int, n_sim: int):
    dt = T/steps
    dBt =  np.column_stack([np.random.normal(0, math.sqrt(dt), n_sim) for _ in range(steps+1)])
    dBt[:, 0] = 0
    Bt = dBt.cumsum(1)

    return Bt[:, :steps], dBt[:, :steps]

def bessel_process(dimensions: int, Z0: float, T: int, steps: int, n_sim: int) -> (np.array, np.array):
    summation = np.zeros([n_sim, steps])
    
    for _ in range(dimensions):
        Bt, dBt = brownian_motion(T, steps, n_sim) 
        summation += Bt * dBt * 2
        
    dRt = summation + dimensions * T/steps
    Rt = dRt.cumsum(1)
    
    return Rt, np.sqrt(Rt)
```

The problem is that the resulting process has a mean 3 times the expected and I can't figure out where is the error in my code.

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.