diff --git a/benchmarks/cupy/examples/finance/__init__.py b/benchmarks/cupy/examples/finance/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/benchmarks/cupy/examples/finance/black_scholes.py b/benchmarks/cupy/examples/finance/black_scholes.py
new file mode 100644
index 0000000..5551cc4
--- /dev/null
+++ b/benchmarks/cupy/examples/finance/black_scholes.py
@@ -0,0 +1,91 @@
+from benchmarks.numpy.common import Benchmark
+from benchmarks.utils import sync
+from benchmarks.utils.helper import parameterize
+
+import cupy
+
+
+black_scholes_kernel = cupy.ElementwiseKernel(
+    'T s, T x, T t, T r, T v',  # Inputs
+    'T call, T put',  # Outputs
+    '''
+    const T sqrt_t = sqrt(t);
+    const T d1 = (log(s / x) + (r + v * v / 2) * t) / (v * sqrt_t);
+    const T d2 = d1 - v * sqrt_t;
+    const T cnd_d1 = get_cumulative_normal_distribution(d1);
+    const T cnd_d2 = get_cumulative_normal_distribution(d2);
+    const T exp_rt = exp(- r * t);
+    call = s * cnd_d1 - x * exp_rt * cnd_d2;
+    put = x * exp_rt * (1 - cnd_d2) - s * (1 - cnd_d1);
+    ''',
+    'black_scholes_kernel',
+    preamble='''
+    __device__
+    inline T get_cumulative_normal_distribution(T x) {
+        const T A1 = 0.31938153;
+        const T A2 = -0.356563782;
+        const T A3 = 1.781477937;
+        const T A4 = -1.821255978;
+        const T A5 = 1.330274429;
+        const T RSQRT2PI = 0.39894228040143267793994605993438;
+        const T W = 0.2316419;
+        const T k = 1 / (1 + W * abs(x));
+        T cnd = RSQRT2PI * exp(- x * x / 2) *
+            (k * (A1 + k * (A2 + k * (A3 + k * (A4 + k * A5)))));
+        if (x > 0) {
+            cnd = 1 - cnd;
+        }
+        return cnd;
+    }
+    ''',
+)
+
+
+@sync
+@parameterize([('n_options', [1, 10000000])])
+class BlackScholes(Benchmark):
+    def setup(self, n_options):
+        def rand_range(m, M):
+            samples = cupy.random.rand(n_options)
+            return (m + (M - m) * samples).astype(cupy.float64)
+        self.stock_price = rand_range(5, 30)
+        self.option_strike = rand_range(1, 100)
+        self.option_years = rand_range(0.25, 10)
+        self.risk_free = 0.02
+        self.volatility = 0.3
+
+    def time_black_scholes(self, n_options):
+        s, x, t = self.stock_price, self.option_strike, self.option_years
+        r, v = self.risk_free, self.volatility
+
+        sqrt_t = cupy.sqrt(t)
+        d1 = (cupy.log(s / x) + (r + v * v / 2) * t) / (v * sqrt_t)
+        d2 = d1 - v * sqrt_t
+
+        def get_cumulative_normal_distribution(x):
+            A1 = 0.31938153
+            A2 = -0.356563782
+            A3 = 1.781477937
+            A4 = -1.821255978
+            A5 = 1.330274429
+            RSQRT2PI = 0.39894228040143267793994605993438
+            W = 0.2316419
+
+            k = 1 / (1 + W * cupy.abs(x))
+            cnd = RSQRT2PI * cupy.exp(-x * x / 2) * (
+                k * (A1 + k * (A2 + k * (A3 + k * (A4 + k * A5)))))
+            cnd = cupy.where(x > 0, 1 - cnd, cnd)
+            return cnd
+
+        cnd_d1 = get_cumulative_normal_distribution(d1)
+        cnd_d2 = get_cumulative_normal_distribution(d2)
+
+        exp_rt = cupy.exp(- r * t)
+        call = s * cnd_d1 - x * exp_rt * cnd_d2
+        put = x * exp_rt * (1 - cnd_d2) - s * (1 - cnd_d1)
+        return call, put
+
+    def time_black_scholes_kernel(self, n_options):
+        black_scholes_kernel(
+            self.stock_price, self.option_strike, self.option_years,
+            self.risk_free, self.volatility)
diff --git a/benchmarks/cupy/examples/finance/monte_carlo.py b/benchmarks/cupy/examples/finance/monte_carlo.py
new file mode 100644
index 0000000..44b98f9
--- /dev/null
+++ b/benchmarks/cupy/examples/finance/monte_carlo.py
@@ -0,0 +1,103 @@
+from benchmarks.numpy.common import Benchmark
+from benchmarks.utils import sync
+from benchmarks.utils.helper import parameterize
+
+import cupy
+
+
+monte_carlo_kernel = cupy.ElementwiseKernel(
+    'T s, T x, T t, T r, T v, int32 n_samples, int32 seed', 'T call',
+    '''
+    // We can use special variables i and _ind to get the index of the thread.
+    // In this case, we used an index as a seed of random sequence.
+    uint64_t rand_state[2];
+    init_state(rand_state, i, seed);
+
+    T call_sum = 0;
+    const T v_by_sqrt_t = v * sqrt(t);
+    const T mu_by_t = (r - v * v / 2) * t;
+
+    // compute the price of the call option with Monte Carlo method
+    for (int i = 0; i < n_samples; ++i) {
+        const T p = sample_normal(rand_state);
+        call_sum += get_call_value(s, x, p, mu_by_t, v_by_sqrt_t);
+    }
+    // convert the future value of the call option to the present value
+    const T discount_factor = exp(- r * t);
+    call = discount_factor * call_sum / n_samples;
+    ''',
+    preamble='''
+    typedef unsigned long long uint64_t;
+
+    __device__
+    inline T get_call_value(T s, T x, T p, T mu_by_t, T v_by_sqrt_t) {
+        const T call_value = s * exp(mu_by_t + v_by_sqrt_t * p) - x;
+        return (call_value > 0) ? call_value : 0;
+    }
+
+    // Initialize state
+    __device__ inline void init_state(uint64_t* a, int i, int seed) {
+        a[0] = i + 1;
+        a[1] = 0x5c721fd808f616b6 + seed;
+    }
+
+    __device__ inline uint64_t xorshift128plus(uint64_t* x) {
+        uint64_t s1 = x[0];
+        uint64_t s0 = x[1];
+        x[0] = s0;
+        s1 = s1 ^ (s1 << 23);
+        s1 = s1 ^ (s1 >> 17);
+        s1 = s1 ^ s0;
+        s1 = s1 ^ (s0 >> 26);
+        x[1] = s1;
+        return s0 + s1;
+    }
+
+    // Draw a sample from an uniform distribution in a range of [0, 1]
+    __device__ inline T sample_uniform(uint64_t* state) {
+        const uint64_t x = xorshift128plus(state);
+        // 18446744073709551615 = 2^64 - 1
+        return T(x) / T(18446744073709551615);
+    }
+
+    // Draw a sample from a normal distribution with N(0, 1)
+    __device__ inline T sample_normal(uint64_t* state) {
+        T x = sample_uniform(state);
+        T s = T(-1.4142135623730950488016887242097);  // = -sqrt(2)
+        if (x > 0.5) {
+            x = 1 - x;
+            s = -s;
+        }
+        T p = x + T(0.5);
+        return s * erfcinv(2 * p);
+    }
+    ''',
+)
+
+
+@sync
+@parameterize([('n_options', [1, 1000]),
+               ('n_samples_per_thread', [1, 1000]),
+               ('n_threads_per_option', [1, 100000])])
+class MonteCarlo(Benchmark):
+    def setup(self, n_options, n_samples_per_thread, n_threads_per_option):
+        def rand_range(m, M):
+            samples = cupy.random.rand(n_options)
+            return (m + (M - m) * samples).astype(cupy.float64)
+        self.stock_price = rand_range(5, 30)
+        self.option_strike = rand_range(1, 100)
+        self.option_years = rand_range(0.25, 10)
+        self.risk_free = 0.02
+        self.volatility = 0.3
+        self.seed = 0
+
+    def time_compute_option_prices(
+            self, n_options, n_samples_per_thread, n_threads_per_option):
+
+        call_prices = cupy.empty(
+            (n_options, n_threads_per_option), dtype=cupy.float64)
+        monte_carlo_kernel(
+            self.stock_price[:, None], self.option_strike[:, None],
+            self.option_years[:, None], self.risk_free, self.volatility,
+            n_samples_per_thread, self.seed, call_prices)
+        return call_prices.mean(axis=1)