Reflected Brownian motion

The marginal distribution (transient distribution) of a one-dimensional Brownian motion starting at 0 restricted to positive values (a single reflecting barrier at 0) with drift μ and variance σ² is

${\displaystyle \mathbb {P} (Z(t)\leq z)=\Phi \left({\frac {z-\mu t}{\sigma t^{1/2}}}\right)-e^{2\mu z/\sigma ^{2}}\Phi \left({\frac {-z-\mu t}{\sigma t^{1/2}}}\right)}$

for all t ≥ 0, (with Φ the cumulative distribution function of the normal distribution) which yields (for μ < 0) when taking t → ∞ an exponential distribution[2]

${\displaystyle \mathbb {P} (Z<z)=1-e^{2\mu z/\sigma ^{2}}.}$

For fixed t, the distribution of Z(t) coincides with the distribution of the running maximum M(t) of the Brownian motion,

${\displaystyle Z(t)\sim M(t)=\sup _{s\in [0,t]}X(s).}$

But be aware that the distributions of the processes as a whole are very different. In particular, M(t) is increasing in t, which is not the case for Z(t).

The heat kernel for reflected Brownian motion at ${\displaystyle p_{b}}$:

${\displaystyle f(x,p_{b})={\frac {e^{-((x-u)/a)^{2}/2}+e^{-((x+u-2p_{b})/a)^{2}/2}}{a(2\pi )^{1/2}}}}$

For the plane above ${\displaystyle x\geq p_{b}}$

The stationary distribution of a reflected Brownian motion in multiple dimensions is tractable analytically when there is a product form stationary distribution,[9] which occurs when the process is stable and[10]

${\displaystyle 2\Sigma =RD+DR'}$

where D = diag(Σ). In this case the probability density function is[7]

${\displaystyle p(z_{1},z_{2},\ldots ,z_{d})=\prod _{k=1}^{d}\eta _{k}e^{-\eta _{k}z_{k}}}$

where η_k = 2μ_kγ_k/Σ_kk and γ = R⁻¹μ. Closed-form expressions for situations where the product form condition does not hold can be computed numerically as described below in the simulation section.

In one dimension the simulated process is the absolute value of a Wiener process. The following MATLAB program creates a sample path.[11]

% rbm.m
n = 10^4; h=10^(-3); t=h.*(0:n); mu=-1;
X = zeros(1, n+1); M=X; B=X;
B(1)=3; X(1)=3;
for k=2:n+1
    Y = sqrt(h) * randn; U = rand(1);
    B(k) = B(k-1) + mu * h - Y;
    M = (Y + sqrt(Y ^ 2 - 2 * h * log(U))) / 2;
    X(k) = max(M-Y, X(k-1) + h * mu - Y);
end
subplot(2, 1, 1)
plot(t, X, 'k-');
subplot(2, 1, 2)
plot(t, X-B, 'k-');

The error involved in discrete simulations has been quantified.[12]

QNET allows simulation of steady state RBMs.[13][14][15]