• C++ Programming for Financial Engineering
    Highly recommended by thousands of MFE students. Covers essential C++ topics with applications to financial engineering. Learn more Join!
    Python for Finance with Intro to Data Science
    Gain practical understanding of Python to read, understand, and write professional Python code for your first day on the job. Learn more Join!
    An Intuition-Based Options Primer for FE
    Ideal for entry level positions interviews and graduate studies, specializing in options trading arbitrage and options valuation models. Learn more Join!

Monte Carlo Simulation to estimate e

Joined
11/5/14
Messages
294
Points
53
Hi all,

Some of you might already know, how to write a Monte-Carlo simulation to estimate [imath]\pi[/imath]. We use the uniform random number generator to produce a pair of pseudo-random numbers [imath](x_i,y_i)[/imath] in [imath][0,1][/imath]. After [imath]N[/imath] draws, we find the proportion of points lying in the unit circle, that is, those with [imath]x_i^2 + y_i^2 \leq 1[/imath] and divide this by the total number of draws. If [imath]N[/imath] is large, this ratio approaches the geometrical probability of the event that a random point drawn from the unit square lies in a unit quarter-circle: [imath]\pi/4[/imath].

While browsing Math SE, I recently came across a somewhat more challenging puzzle/tasks. How do you estimate the Euler's constant [imath]e \approx 2.71828\ldots [/imath] using the same generator?

My initial idea : If [imath]X[/imath] is [imath]Expo(1)[/imath] random variable, its CDF is [imath]F_X(x) = 1 - e^{-x}[/imath]. By the inverse transform method, [imath]F_X^{-1}(U)[/imath] is an exponential random variable. I generated a sample of [imath]U[0,1][/imath] random numbers, and applied the inverse transform. Since, [imath]F_X(1)=1-e^{-1}[/imath], we can find a count of the number of sample points less than or equal to unity. Their proportion should approach [imath]1-e^{-1}[/imath].

C++:
import numpy as np
import math

def generate_e(N):
    sum = 0

    # Draw N samples from a uniform distribution
    U = np.random.uniform(low=1e-9,high=1.0,size=N)

    # Using the inverse transform technique, log(1/(\lambda(1-u)))~Expo(\lambda) distribution
    lambda_ = 1.0
    X = -np.log(lambda_ *(1.0 - U))

    X2 = X[X<=1.0]
    sum = len(X2)

    e = 1 / (1 - (sum / N))
    return e

estimate = generate_e(10**6)
print(estimate)

However, this approach is really bad (very slow convergence). Also, you aren't allowed to use log and circular functions. Here's the Stats SE thread, that discusses a much simpler and cool idea. (Also, I still need to prove that the random variable [imath]\xi[/imath] defined in the first answer to the OP's question has expectation [imath]e[/imath]. It'll make for an interesting exercise.)
 
Last edited:
V2 u(x) = exp(x) satisfies ODE

du/dx = u
u(0) = 1

Use FDM and compute up to x = 1. QED
(Padé approximant, essentially)

I looked up Pade's approximant on Wiki; so [imath]P_N^M(x)[/imath] is essentially a quotient of two polynomials of degree [imath]M,N[/imath] whose power series ties up upto the first [imath]M+N+1[/imath] terms with the power series of the function.

For the function [imath]f(x)=e^x[/imath], I computed

[math] P_3^3(x) = \frac{1+\frac{1}{2}x + \frac{1}{10}x^2 + \frac{1}{120}x^3}{1-\frac{1}{2}x + \frac{1}{10}x^2 - \frac{1}{120}x^3} [/math]
I get [imath]P_3^3(1) = 2.71830986[/imath], having a relative error [imath]1.03 \times 10^{-5}[/imath], still better than Taylor's series truncated upto terms of degree [imath]6[/imath]; a relative error [imath]8.32 \times 10^{-5}[/imath].

It's like the Pade's approximant of order [imath]3[/imath] is trying to kill [imath]6[/imath] terms of the Taylor series. Pretty impressive speedup!

Btw Daniel, the original puzzle was to use a non-deterministic way to approximate [imath]e[/imath].
 
Last edited:
I looked up Pade's approximant on Wiki; so [imath]P_N^M(x)[/imath] is essentially a quotient of two polynomials of degree [imath]M,N[/imath] whose power series ties up upto the first [imath]M+N+1[/imath] terms with the power series of the function.

For the function [imath]f(x)=e^x[/imath], I computed

[math] P_3^3(x) = \frac{1+\frac{1}{2}x + \frac{1}{10}x^2 + \frac{1}{120}x^3}{1-\frac{1}{2}x + \frac{1}{10}x^2 - \frac{1}{120}x^3} [/math]
I get [imath]P_3^3(1) = 2.71830986[/imath], having a relative error [imath]1.03 \times 10^{-5}[/imath], still better than Taylor's series truncated upto terms of degree [imath]6[/imath]; a relative error [imath]8.32 \times 10^{-5}[/imath].

It's like the Pade's approximant of order [imath]3[/imath] is trying to kill [imath]6[/imath] terms of the Taylor series. Pretty impressive speedup!

Btw Daniel, the original puzzle was to use a non-deterministic way to approximate [imath]e[/imath].
That is Pade(q, p), more generally.
Very important for PDE and ODE

Pade(0,1) implicit Euler
Pade(1,0) explicit Euler
Pade(1,1) Crank Nicolson

Gives a great foothold into FDM. You could run them as well to investigate accuracy O(dt^(p+q+1)) as dt -> 0.
 
Last edited:
That is Pade(q, p), more generally.
Very important for PDE and ODE

Pade(0,1) implicit Euler
Pade(1,0) explicit Euler
Pade(1,1) Crank Nicolson

Gives a great foothold into FDM. You could run them as well to investigate accuracy O(dt^(p+q+1)) as dt -> 0.
Simple foward Euler and backward Euler

C++:
import numpy as np
import math
import matplotlib.pyplot as plt

plt.style.use('seaborn-poster')

f = lambda t, u: u          # ODE
h = 0.10                    # Step-size
A = 0.0
B = 1.0                     # Region
u_0 = 1                     # Initial condition

def forwardEuler(f,A,B,h,u_0):
    t = np.arange(A, B + h, h)  # Define the mesh points
    u = np.zeros(len(t))
    u[0] = u_0

    for k in range(0,len(t)-1):
        u[k+1] = u[k] + h*f(t[k],u[k])

    return u,t

def backwardEuler(f,A,B,h,u_0):
    t = np.arange(A, B + h, h)  # Define the mesh points
    u = np.zeros(len(t))
    u[0] = u_0

    for k in range(0,len(t)-1):
        u[k+1] = u[k]/(1-h)

    return u,t

The error is clearly [imath]O(h)[/imath].

Absolute error in fwd difference for h = 0.1, err = 0.124539368359045
Absolute error in fwd difference for h = 0.05, err = 0.06498412331462378
Absolute error in bkwd difference for h = 0.1, err = 0.14969016233339527
Absolute error in bkwd difference for h = 0.05, err = 0.07122798905721561
 

Attachments

  • ODE_solver.png
    ODE_solver.png
    104.8 KB · Views: 8
Very nice.
The Crank Nicolcon Pade(1,1) is 2nd order and is used a _lot_ in finance.

Now, we can discretise Black Scholes pde(S,t) in S to get an ODE(t) and use the above ODE solvers and others (Boost C++ odeint, Python odeitnt).
aka Method of Lines (MOL).
 

Attachments

  • Dissertation_1381380.pdf
    850 KB · Views: 14
Back
Top