Examples#
Every example below is a complete, runnable program. The code shown is read directly from the same files the playground loads, so what you read here is exactly what runs — click Open in playground under any example to try it.
The first four are period/iteration models, where a variable is recomputed each iteration. The last three are event models, driven by a scheduled clock — including the two written as agents, which is a way of expressing an event model as one participant’s story rather than as a set of handlers.
Estimating Pi#
The classic Monte Carlo estimate: throw random points at a square and count
how many land inside the inscribed circle. The ratio approaches pi/4.
Every pi variable is recomputed once per iteration: two coordinates are
drawn, circle_points uses a guarded expression to increment only on a
hit, and total_points counts every throw. Because est is declared
obs_r, each of the 50 runs contributes one observation, and the boxplot
shows how much a single estimate can vary.
# Estimating Pi
#
# Throw random points at a square and count how many land inside the
# inscribed circle. The ratio of hits to throws approaches pi/4, so four
# times that ratio estimates pi.
# model
pi x := :randu(-1.0, 1.0) # a new point every iteration
pi y := :randu(-1.0, 1.0)
pi circle_points :=
| x^2 + y^2 <= 1.0 -> circle_points + 1
| _ -> circle_points
pi total_points := total_points + 1
obs_r est := 4.0 * (circle_points / total_points)
exec
# execution
P := 500_000 # iterations in one run
R := 50 # number of runs
:randseed(1043)
run() {P, R}
# results
:printf("P: [%d], R: [%d]\n", P, R)
:print(:boxplot(est)) # R observations, one estimate per run
:print("median estimate of Pi", :median(est))
Geometric Brownian Motion#
The standard model of a stock price: a constant drift plus a random shock each day, compounded.
This is the simplest form of a recurrence — the [0] = line seeds the
variable, and every iteration after that refers to the variable’s own
previous value. Observing price with obs_r captures where each of the
10,000 paths ended up after a year of trading days.
# Geometric Brownian Motion
#
# The standard model of a stock price: constant drift plus random shocks,
# compounded daily. Unlike idx500.eb there are no cashflows -- the whole
# path is driven by the return process itself.
# model
S0 := 100.0 # initial price, dollars
mu := 0.05 # expected annual return (drift)
sigma := 0.20 # annual volatility
dt := 1.0 / 252.0 # one trading day, in years
pi price: dbl =
[0] = S0
price * :exp((mu - 0.5*sigma^2)*dt + sigma*:sqrt(dt)*:randn(0.0, 1.0))
obs_r price
exec
# execution
P := 252 # trading days in one run (one year)
R := 10_000 # number of runs
:randseed(1043)
run() {P, R}
# results
:printf("P: [%d], R: [%d]\n", P, R)
:print(:boxplot(price)) # R observations, the final price of each run
:printf("median final price: %0.2f\n", :median(price))
Index 500 Portfolio#
The same compounding idea put to work on a real question: if you invest $10,000 and add another $10,000 every year for 30 years, how likely are you to reach $1,000,000?
The recurrence now carries a cashflow as well as a return. More importantly,
this is the first example with scenario blocks: each one overrides model
parameters, and run takes both the scenario and a label, so the two
boxplots come out identified. Comparing a calm market against a volatile one
costs five lines rather than a second copy of the model.
# Index 500 Portfolio
#
# Determine the probability of achieving a financial goal by investing in an
# Index 500 fund, adding money at the end of each period. Extends a demo by
# Matt Macarty, "Basic Monte Carlo Simulation of a Stock Portfolio || Python
# Programming": https://www.youtube.com/watch?v=A0J0VAHzIxc
# model
rate := 0.07 # expected return
sd := 0.07 # expected standard deviation
goal := 1_000_000 # desired outcome, dollars
initial := 10_000 # initial investment, dollars
added := 10_000 # added investment at end of each year, dollars
pi mkt_return := :randn(rate, sd)
pi value: dbl =
[0] = initial
value * (1 + mkt_return) + added
pi reached_goal := | value >= goal -> true
| _ -> false
obs_r value
obs_r reached_goal
exec
# execution
scenario base
# use initial values set in model
end
scenario high_sd
rate = 0.09
sd = 0.17
end
P := 30 # periods (years) in one run
R := 5_000 # number of runs
:randseed(1043)
# results
run(base, "default volatility") {P, R}
:print(:boxplot(value)) # R observations, the final value of each run
:print("median result", :median(value))
:printf("reached goal: %0.2f%%\n", 100 * :probability(reached_goal))
run(high_sd, "high volatility") {P, R}
:print(:boxplot(value))
:print("median result", :median(value))
:printf("reached goal: %0.2f%%\n", 100 * :probability(reached_goal))
Population Extinction#
A branching process: each individual independently leaves 0, 1, or 2 offspring. Mean offspring here is 1.2, so the population grows on average — yet it still dies out about a third of the time, and 1/3 is exactly the analytic answer.
Three new things appear. A helper fn does the per-generation tally.
stop_r ends a run early once the population hits zero, since nothing can
happen afterwards. And extinction_gen reads :p inside the model to
stamp which generation extinction happened in — legitimate use of the
counter, because it happens while the loop is running.
# Population Extinction
#
# Estimate the probability of extinction in a branching process, where each
# individual independently leaves 0, 1, or 2 offspring. For the distribution
# below the analytic answer is 1/3. See the Wikipedia entry on "Branching
# Processes".
# model
# Offspring distribution, as cumulative cutoffs on :randu(0, 1):
# 0 children with p = 0.1, 1 child with p = 0.6, 2 children with p = 0.3
# Mean offspring is 1.2 -- above 1, so the population can grow, yet a third
# of runs still die out.
p1 := 0.1
p2 := 0.7
fn sum(p: int, parm1: dbl, parm2: dbl): int # offspring of one generation
total := 0
for _ = 1, p do
r := :randu(0, 1)
k := | r < parm1 -> 0
| r < parm2 -> 1
| _ -> 2
total = total + k
end
return total
end
pi pop :=
[0] = 1
| pop == 0 -> 0
| _ -> sum(pop, p1, p2)
pi extinct := | pop == 0 -> true
| _ -> false
pi extinction_gen := | extinct -> :p
| _ -> nil
obs_r extinct
obs_r extinction_gen
stop_r(extinct == true, "population extinct") # nothing can happen after
exec
# execution
P := 25 # generations in one run
R := 500_000 # number of runs
:randseed(1043)
run() {P, R}
# results
:printf("extinction probability: %0.2f%%\n", 100 * :probability(extinct))
:printf("median extinction generation: %0.2f\n", :median(extinction_gen))
Gambler’s Ruin#
A gambler bets $1 at a time, starting with $50 and quitting at $0 or $100. With a fair coin the probability of ruin is 0.5; shading the odds slightly against them, to 0.45, makes ruin nearly certain.
This model is driven by events rather than iterations. state variables
persist across a run and reset at the start of each one, an eh handler
reschedules itself with sched,
and the run ends naturally when the handler stops rescheduling and the event
queue empties — no explicit halt. The file’s closing notes work through why
this is a state machine rather than a time-based model.
# Gambler's Ruin
#
# A gambler bets $1 at a time, starting at $50 and quitting at $0 or $100.
# With a fair coin the probability of ruin is 0.50; shading the odds to 0.45
# makes ruin all but certain. Each run ends naturally when the event queue
# empties, which happens once wealth reaches either boundary.
# model
start := 50
target := 100
p := 0.5
state wealth: int = start
state ruined: int = 0
eh step()
if :randu() < p then
wealth = wealth + 1
else
wealth = wealth - 1
end
if wealth > 0 and wealth < target then
sched(1, step)
end
if wealth <= 0 then
ruined = 1
end
end
sched(0, step) # seed the chain; fires once at t=0 of each run
obs_r ruined
exec
# execution
scenario fair
# use initial values set in model
end
scenario unfair
p = 0.45
end
T := 100_000 # time limit for one run
R := 1_000 # number of runs
:randseed(1043)
# results
run(fair, "fair") {T, R} # p(ruin) ~ 0.50
:print(:mean(ruined))
run(unfair, "unfair") {T, R} # p(ruin) ~ 0.9999
:print(:mean(ruined))
# A few design notes:
#
# Natural termination. When wealth hits 0 or target, the handler doesn't
# reschedule itself. The event queue empties and the run ends -- no halt or
# stop_r needed.
#
# sched(0, step) at module level. This is the idiomatic way to seed an EH
# simulation: fires once at t=0 at the start of each run, starting the chain.
#
# Why this is "not time-based." The handler does use sched(1, step) so :t
# counts steps, but :t is never observed. The outcome of interest is ruined --
# a state flag, not a time. It's a state machine {playing -> playing, playing
# -> ruined, playing -> success} driven by events, and what we measure is
# which absorbing state each run ends in.
#
# Theoretical values. For p=0.5, P(ruin) = start/target = 0.5. For p=0.45,
# P(ruin) = r^50/(1 + r^50) where r = q/p = 0.55/0.45 ~ 1.22, giving
# r^50 ~ 22900 and P(ruin) ~ 0.9999.
Bank Teller: Two Customers#
Two customers arrive at a bank with a single teller. The second has to wait for the first to finish.
An ag body reads as the story of one customer from arrival to departure,
and the simulation is what happens when several such stories overlap. This is
the same event machinery as Gambler’s Ruin above, written from the
participant’s point of view instead of the handler’s. The teller is a
resource (rc), and claim waits until it is free — which is the entire
source of the queueing behavior. init stages the arrivals.
# Bank Teller: Two Customers
#
# Two customers arrive at a bank with a single teller. The second has to
# wait for the first to finish -- that wait is the whole of the queueing
# behavior, and it comes entirely from claim().
# model
rc teller := 1 # resource; number of tellers
state last_leave: dbl = 0.0 # time the last customer left
ag customer(name: str, res: rsc)
:printf("%s arrives at t=%0.2f\n", name, :t)
claim(res)
:printf("%s starts at t=%0.2f\n", name, :t)
timeout(4)
end
:printf("%s leaves at t=%0.2f\n", name, :t)
last_leave = :max(last_leave, :t)
end
obs_r last_leave
init
customer("Customer 1", teller)
timeout(1)
customer("Customer 2", teller)
end
exec
# execution
T := 100 # time limit for one run
R := 1 # number of runs
run() {T, R}
# results
:print(:mean(last_leave))
Bank Teller: N Customers#
The same bank, opened up: N customers arrive a minute apart, and service
times are now random rather than a flat four minutes.
Nothing about the agent changed — only how many are created. The scenario
pair raises N from 20 to 100 against the same single teller, which is
enough to push the queue from manageable to hopeless.
# Bank Teller: N Customers
#
# The same bank opened up: N customers arrive a minute apart, and service
# times are now random rather than a flat four minutes. Nothing about the
# agent changes -- only how many of them there are.
# model
N := 20 # number of customers
rc teller := 1 # resource; number of tellers
state last_leave: dbl = 0.0 # time the last customer left
ag customer(id: int, res: rsc)
:printf("Customer %d arrives at t=%0.2f\n", id, :t)
claim(res)
:printf("Customer %d starts at t=%0.2f\n", id, :t)
timeout(:rande(4))
end
:printf("Customer %d leaves at t=%0.2f\n", id, :t)
last_leave = :max(last_leave, :t)
end
obs_r last_leave
init
for i = 1, N do
customer(i, teller)
timeout(1)
end
end
exec
# execution
scenario base
# use initial values set in model
end
scenario high_volume
N = 100
end
T := 100 # time limit for one run
R := 1 # number of runs
:randseed(1043)
# results
run(base, "base") {T, R}
:print(:mean(last_leave))
run(high_volume, "high volume") {T, R}
:print(:mean(last_leave))