Probability of a 2D Brownian motion exiting a circle of radius r before time t.

Question

Suppose P(x,y)_t is following 2D Brownian motion. That is x,y are IID N(0,t). If P_0 = (0,0) then I want to find the probability distribution of T = inf{ t : Norm(P(x,y)_t) > 1}. First step was that I did the simulations to find the empirical answer. Then, I thought I need to find P(X^2 +Y^2 > 1) which should be P(ChiSq_df_2 > 1) = 0.6 which wasn't the answer coming through simulations (0.88)(Time = 1).

After spending some time on the internet I found out that the radial length of a 2d Brownian motion follows a stochastic differential equation called Bessel Process.

The Squared Bessel Process is this : dX_t = 2(X_t)^(1/2) dB_t + 2dt.

The good thing was simulated squared Bessel process and the actual simulation agreed with each other.

But, I also read that the squared Bessel process follows Non Central Chi Squared distribution, which again brings me to where I started.

Please help me find the distribution of X_t which should be chi square but simulations don't agree.

I am attaching my simulation codes :

This one is for simulating the 2d Brownian motion

import matplotlib.pyplot as plt

def simulate_brownian_exit_time(r, dt=0.01, num_simulations=1000):
    exit_times = []

    for _ in range(num_simulations):
        x, y = 0.0, 0.0  # Start at the origin
        time = 0

        while np.sqrt(x**2 + y**2) < r:
            dx, dy = np.random.normal(0, np.sqrt(dt), 2)  # 2D Brownian step
            x += dx
            y += dy
            time += dt

        exit_times.append(time)

    return exit_times

def plot_exit_time_distribution(exit_times):
    plt.hist(exit_times, bins=30, color='skyblue', edgecolor='black')
    plt.title('Distribution of Exit Times for 2D Brownian Motion')
    plt.xlabel('Exit Time')
    plt.ylabel('Frequency')
    plt.show()

# Parameters
r = 1  # Circle radius
dt = 0.01  # Time step
num_simulations = 10000  # Number of simulations

# Simulate and plot
exit_times = simulate_brownian_exit_time(r, dt, num_simulations)
t = 1
exits_before_t = len([time for time in exit_times if time < t])
probability_before_t = exits_before_t / num_simulations
print(f"Probability of exit before time {t}: {probability_before_t}")

average_exit_time = np.mean(exit_times)
print(f"Average exit time: {average_exit_time} units")
plot_exit_time_distribution(exit_times)

and this is the code for simulating squared Bessel process :

import matplotlib.pyplot as plt

# Parameters
dt = 0.01  # Time step
X0 = 0.0  # Initial value of the squared Bessel process
r_squared = 1  # Target squared radius to compare against
specified_time_t = 1  # Specific time t to evaluate probability of exiting before

# Monte Carlo Parameters
simulations = 10000  # Number of Monte Carlo simulations
exit_times = []  # List to store exit times

# Simulation with progress update for average and probability
for sim in range(simulations):
    X = X0
    time = 0
    while X <= r_squared:
        Z = np.random.normal(0, 1)  # Standard normal random variable
        X += 2 * np.sqrt(abs(X)) * np.sqrt(dt) * Z + 2 * dt
        time += dt
    exit_times.append(time)  # Store the exit time for this simulation

    # Progress update after every 10% of simulations
    if (sim + 1) % (simulations // 10) == 0:
        # Calculate current average exit time
        current_average_exit_time = np.mean(exit_times)
        # Calculate current probability of exiting before specified time t
        current_exits_before_t_count = np.sum(np.array(exit_times) < specified_time_t)
        current_probability_before_t = current_exits_before_t_count / (sim + 1)
        print(f"After {(sim + 1)} simulations, current average exit time: {current_average_exit_time:.4f}, current probability before t={specified_time_t}: {current_probability_before_t:.4f}")

# Final calculations
average_exit_time = np.mean(exit_times)
probability_before_t = np.sum(np.array(exit_times) < specified_time_t) / simulations

print(f"\nFinal Average exit time: {average_exit_time:.4f} units")
print(f"Final Probability of exit before time t={specified_time_t}: {probability_before_t:.4f}")

# Plotting the exit time distribution
plt.hist(exit_times, bins=50, color='skyblue', edgecolor='black')
plt.title('Distribution of Exit Times for Squared Bessel Process')
plt.xlabel('Exit Time')
plt.ylabel('Frequency')
plt.show()
```