0Pricing
R Academy · Lesson

Monte Carlo Simulation Basics

Estimate pi and other quantities by repeated random sampling.

Monte Carlo Simulation Basics is a free R Academy lesson on CoddyKit — lesson 3 of 4. You can read the complete lesson below for free — then practise it hands-on in the browser with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, one of 4 lessons in the course, and your progress syncs across the web and the CoddyKit app.

What Is Monte Carlo?

Monte Carlo simulation uses repeated random sampling to estimate quantities that are hard to compute analytically. Named after the casino in Monaco, it's used in finance, physics, statistics, and AI for integration, optimization, and inference.

# Core idea: approximate deterministic quantities
# using random sampling and the Law of Large Numbers

# Example: estimate probability that sum of two dice > 8
set.seed(42)
n <- 10000
die1 <- sample(1:6, n, replace = TRUE)
die2 <- sample(1:6, n, replace = TRUE)
total <- die1 + die2

# Monte Carlo estimate
mc_estimate <- mean(total > 8)
cat('MC estimate P(sum > 8):', mc_estimate, '\n')

# Exact probability
exact <- sum(outer(1:6, 1:6, '+') > 8) / 36
cat('Exact probability:', exact)

Estimating Pi with Random Points

The classic Monte Carlo demo: randomly place points in a unit square. The fraction falling inside the unit circle estimates π/4. As n → ∞, the estimate converges to π.

set.seed(42)
n <- 100000

# Random points in [-1, 1] x [-1, 1] unit square
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)

# Check if inside unit circle: x^2 + y^2 <= 1
inside <- (x^2 + y^2) <= 1

# Pi estimate: fraction inside * area of square
pi_estimate <- 4 * mean(inside)
cat('Pi estimate:', pi_estimate, '\n')
cat('True pi:', pi, '\n')
cat('Error:', abs(pi_estimate - pi))

Estimating Expected Value

Monte Carlo integration: E[f(X)] ≈ (1/n) Σ f(Xᵢ) where X₁,...,Xₙ are iid draws. Integrating g(x) from a to b: sample x from Uniform(a,b), E[g(X)] ≈ (b-a) * mean(g(samples)).

set.seed(42)
# Estimate integral of sin(x) from 0 to pi
# Exact value = 2
n <- 50000
x_samples <- runif(n, min = 0, max = pi)
integrand_vals <- sin(x_samples)

# MC estimate = (b-a) * mean(f(x))
mc_integral <- pi * mean(integrand_vals)  # (pi - 0) * mean
cat('MC estimate of integral:', mc_integral, '\n')
cat('Exact value: 2\n')
cat('Error:', abs(mc_integral - 2), '\n')

# Standard error of the estimate
se <- pi * sd(integrand_vals) / sqrt(n)
cat('Standard error:', se)

replicate() for Repeated Simulations

replicate(n, expr) is R's built-in function for running a simulation n times and collecting results. It's cleaner than a for loop for simple simulations and returns a vector or matrix.

set.seed(42)
# Simulate sample mean of 30 N(0,1) draws
# Repeat 10000 times to study sampling distribution
sample_means <- replicate(10000, {
  x <- rnorm(30)  # sample of 30
  mean(x)          # compute mean
})

# Central Limit Theorem: sample mean ~ N(0, 1/sqrt(30))
mean(sample_means)  # ~0
sd(sample_means)    # ~1/sqrt(30) = 0.183

# 95% CI width
diff(quantile(sample_means, c(0.025, 0.975)))

# Compare to theoretical
2 * 1.96 / sqrt(30)

Law of Large Numbers Demo

The Law of Large Numbers says the sample mean converges to the true mean as n increases. Watching this convergence in R demonstrates why Monte Carlo works and how fast estimates stabilize.

set.seed(42)
# Rolling a fair die: true mean = 3.5
n_max <- 10000
rolls <- sample(1:6, n_max, replace = TRUE)
cumulative_means <- cumsum(rolls) / seq_along(rolls)

# Show convergence at different sample sizes
ns <- c(10, 100, 1000, 5000, 10000)
results <- data.frame(
  n = ns,
  mean = cumulative_means[ns],
  error = abs(cumulative_means[ns] - 3.5)
)
print(results)
# Error decreases as n increases

Variance Reduction: Antithetic Variates

Use uniform samples u and their complements (1-u) as an antithetic pair. Their negative correlation reduces estimator variance by up to 50% with the same number of function evaluations.

set.seed(42)
n <- 1000

# Standard MC: estimate E[exp(U)] where U~Uniform(0,1)
# True value = e - 1 = 1.718282
u <- runif(n)
mc_std <- mean(exp(u))

# Antithetic variates: use u AND 1-u
u_anti <- runif(n/2)
mc_anti <- mean((exp(u_anti) + exp(1 - u_anti)) / 2)

cat('True value:', exp(1) - 1, '\n')
cat('Standard MC:', mc_std, '\n')
cat('Antithetic MC:', mc_anti, '\n')

# Variance comparison
var_std  <- var(exp(runif(10000)))
var_anti <- var((exp(runif(5000)) + exp(1 - runif(5000)))/2)
cat('Variance ratio (anti/std):', var_anti/var_std)

Control Variates

A control variate is a function with known expectation, correlated with the target. Subtract it (scaled) to reduce variance. Classic: estimate E[f(X)] using correlated g(X) whose mean is known.

set.seed(42)
n <- 5000

# Estimate E[exp(U)] where U~Uniform(0,1)
# True: e - 1 = 1.71828

# Control variate: g(U) = U, E[U] = 0.5
u <- runif(n)
f_vals <- exp(u)  # target
g_vals <- u       # control variate

# Optimal coefficient c = -Cov(f,g)/Var(g)
c_star <- -cov(f_vals, g_vals) / var(g_vals)

# Control variate estimator
mc_cv <- mean(f_vals + c_star * (g_vals - 0.5))

cat('Standard MC:', mean(f_vals), '\n')
cat('Control variate MC:', mc_cv, '\n')
cat('True value:', exp(1) - 1)

# Variance reduction factor
var(f_vals) / var(f_vals + c_star * (g_vals - 0.5))

Monte Carlo for Probability Estimation

Monte Carlo excels at estimating complex probabilities that are intractable analytically. Simulate the random process many times and compute the fraction of times the event occurs.

set.seed(123)
# Birthday problem: P(at least 2 people share birthday)
# in a group of n people

birtday_collision <- function(n_people) {
  birthdays <- sample(1:365, n_people, replace = TRUE)
  length(birthdays) != length(unique(birthdays))
}

# Estimate for groups of size 10, 23, 50
sizes <- c(10, 23, 50)
for (sz in sizes) {
  p <- mean(replicate(5000, birt_day_collision <- {
    bd <- sample(1:365, sz, replace = TRUE)
    length(bd) != length(unique(bd))
  }))
  cat('n =', sz, ': P(collision) ~', round(p, 3), '\n')
}

Simulating Geometric Brownian Motion

Stock prices are often modeled as Geometric Brownian Motion: S(t+dt) = S(t) * exp((μ - σ²/2)dt + σ√dt * Z), where Z~N(0,1). Monte Carlo generates price paths for option pricing.

set.seed(42)
S0 <- 100     # initial price
mu <- 0.05    # annual drift
sigma <- 0.2  # annual volatility
T <- 1        # 1 year
n_steps <- 252 # daily steps
dt <- T / n_steps

# Simulate one price path
Z <- rnorm(n_steps)
log_returns <- (mu - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z
price_path <- S0 * exp(cumsum(log_returns))

cat('Final price:', round(price_path[n_steps], 2), '\n')
cat('Min price:', round(min(price_path), 2), '\n')
cat('Max price:', round(max(price_path), 2))

MC Option Pricing

Price a European call option using Monte Carlo: simulate many final stock prices, compute payoffs max(S_T - K, 0), then discount the average payoff by e^(-rT).

set.seed(42)
S0 <- 100; K <- 105; r <- 0.05; sigma <- 0.2; T <- 1
n_sim <- 50000

# Final stock prices under risk-neutral measure
Z <- rnorm(n_sim)
ST <- S0 * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)

# Call option payoff
payoff <- pmax(ST - K, 0)

# Discounted expected payoff
call_price <- exp(-r * T) * mean(payoff)
se <- exp(-r * T) * sd(payoff) / sqrt(n_sim)
cat('Call price:', round(call_price, 4), '\n')
cat('95% CI: [', round(call_price - 1.96*se, 4),
    ',', round(call_price + 1.96*se, 4), ']')

Convergence Rate of MC

Monte Carlo converges at rate O(1/√n): halving the error requires 4× more samples. The standard error of the MC estimate is σ/√n where σ is the standard deviation of the integrand.

set.seed(42)
# Demonstrate MC convergence for pi estimation
estimate_pi <- function(n) {
  x <- runif(n, -1, 1)
  y <- runif(n, -1, 1)
  4 * mean(x^2 + y^2 <= 1)
}

# Sample sizes: powers of 10
ns <- 10^(1:5)
estimates <- sapply(ns, function(n) {
  set.seed(42)
  estimate_pi(n)
})

data.frame(
  n = ns,
  pi_estimate = round(estimates, 5),
  error = round(abs(estimates - pi), 5)
)

Quick Check

Test your understanding of Monte Carlo simulation fundamentals.

Recap: Monte Carlo Basics

Key takeaways: Monte Carlo estimates quantities by averaging over many random samples. Use replicate() for repeated simulations. The Law of Large Numbers guarantees convergence. Error scales as 1/√n — quadruple samples to halve error. Variance reduction techniques (antithetic variates, control variates) improve efficiency without more samples. Applications: integration, probability estimation, option pricing, simulation.

set.seed(42)
# Monte Carlo template:
# 1. Define simulation function
simulate_once <- function() {
  x <- runif(1, -1, 1)
  y <- runif(1, -1, 1)
  (x^2 + y^2) <= 1
}

# 2. Replicate many times
n <- 10000
results <- replicate(n, simulate_once())

# 3. Estimate quantity of interest
pi_mc <- 4 * mean(results)

# 4. Quantify uncertainty
se <- 4 * sd(results) / sqrt(n)
cat('Pi:', pi_mc, '+/-', round(1.96*se, 4))

Frequently asked questions

Is the “Monte Carlo Simulation Basics” lesson free?

Yes — the full text of “Monte Carlo Simulation Basics” is free to read here on the web, and the R Academy course includes 4 lessons in total. To practise it interactively (a built-in code editor and a 24/7 AI tutor) and unlock the rest of the R Academy course, upgrade to CoddyKit PRO.

What will I learn in “Monte Carlo Simulation Basics”?

Estimate pi and other quantities by repeated random sampling. You practise R Academy with hands-on code you run directly in the browser, and a 24/7 AI tutor answers your questions as you work through the lesson.

Do I need any experience to start R Academy?

No prior experience is required. R Academy on CoddyKit is structured for beginners through advanced learners; this is — lesson 3 of 4, so you can start here or from the beginning and move at your own pace.

How long does the “Monte Carlo Simulation Basics” lesson take?

Most CoddyKit lessons take about 5–10 minutes. Each one is bite-sized and interactive, so you make steady progress and pick up exactly where you left off across the web and the app.

Can I write and run code in this R Academy lesson?

Yes. Every R Academy lesson includes a built-in code editor, so you write and run real code right in your browser and get instant AI feedback — no local setup required.

All lessons in this course

  1. set.seed() and Reproducibility
  2. Generating Random Distributions
  3. Monte Carlo Simulation Basics
  4. Bootstrap Resampling in R
← Back to R Academy