Why a Markov Chain Is a (Surprisingly) Good Bet for Tomorrow's Weather¶
Last time, we showed that predicting the exact rainfall on a specific day six months out is essentially impossible — the atmosphere is chaotic, and tiny uncertainties snowball within about two weeks. So here's a much smaller question: if it's raining right now, what's the chance it's still raining tomorrow?
For that question, meteorologists (and statisticians modeling all sorts of systems) reach for one of the simplest tools available: a Markov chain, which assumes tomorrow only cares about today — not the whole history that led here. That sounds almost too lazy to work. Let's actually test it, using four years of real daily weather data.
What's a Markov chain?¶
A Markov chain is a system that hops between a fixed set of states, where the probability of the next state depends only on the current state — not on how it got there. This property has a name: the Markov property, or "memorylessness."
Think of a board game token: to decide where it moves next, you only need to know which square it's currently on and roll a die. You don't need to know the path it took to get there. Weather, modeled this way, has just two squares: Wet and Dry.
It's a bold simplification — real weather obviously does have memory (a slow-moving storm system doesn't reset every 24 hours). The question is whether that memory matters enough to bother modeling, or whether "today" alone already captures most of what we need to know. Let's find out.
stateDiagram-v2
[*] --> Dry
Dry --> Dry: stays dry
Dry --> Wet: turns wet
Wet --> Wet: stays wet
Wet --> Dry: turns dry
Two states, four possible transitions -- that's the whole model. No history, no memory of last week's forecast, just "given where I am now, where do I go next." Let's replace those generic arrows with real probabilities.
import re
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from IPython.display import HTML, display
def show(fig, div_id, height=420):
"""Render a self-contained, interactive Plotly figure (hover, zoom, pan)."""
fig.update_layout(
paper_bgcolor="white",
plot_bgcolor="white",
margin=dict(t=60, b=40, l=50, r=20),
)
html = fig.to_html(
full_html=False,
include_plotlyjs="cdn",
default_height=f"{height}px",
div_id=div_id,
)
# Strip the auto-generated SRI hash on the CDN <script> tag -- it's a long
# base64 string that trips secret scanners as a false positive.
html = re.sub(r'\s+integrity="[^"]*"\s+crossorigin="anonymous"', "", html)
display(HTML(html))
The data¶
We'll use four full years (2012–2015) of daily weather observations for Seattle, Washington — 1,461 days of precipitation, temperature, and wind, published as part of the open-source vega-datasets collection (originally sourced from NOAA). We only need one column: daily precipitation, which we'll collapse into a binary Wet (precipitation > 0mm) / Dry (no precipitation) label for each day.
df = pd.read_csv(
"https://raw.githubusercontent.com/vega/vega-datasets/main/data/seattle-weather.csv"
)
df["date"] = pd.to_datetime(df["date"])
df = df.sort_values("date").reset_index(drop=True)
df["wet"] = (df["precipitation"] > 0).astype(int) # 1 = wet day, 0 = dry day
STATE_NAMES = ["Dry", "Wet"]
df[["date", "precipitation", "wet"]].head()
| date | precipitation | wet | |
|---|---|---|---|
| 0 | 2012-01-01 | 0.0 | 0 |
| 1 | 2012-01-02 | 10.9 | 1 |
| 2 | 2012-01-03 | 0.8 | 1 |
| 3 | 2012-01-04 | 20.3 | 1 |
| 4 | 2012-01-05 | 1.3 | 1 |
Building the chain from real transitions¶
For every pair of consecutive days, we ask: what state was it today, and what state did it become tomorrow? Counting up all 1,460 of those transitions and normalizing each row gives us the chain's transition matrix — the probability of moving from any state to any other.
states = df["wet"].values
today, tomorrow = states[:-1], states[1:]
counts = np.zeros((2, 2))
for s0, s1 in zip(today, tomorrow):
counts[s0, s1] += 1
trans_prob = counts / counts.sum(axis=1, keepdims=True)
pd.DataFrame(
trans_prob,
index=[f"today: {s}" for s in STATE_NAMES],
columns=[f"tomorrow: {s}" for s in STATE_NAMES],
).round(3)
| tomorrow: Dry | tomorrow: Wet | |
|---|---|---|
| today: Dry | 0.756 | 0.244 |
| today: Wet | 0.327 | 0.673 |
fig = go.Figure(
data=go.Heatmap(
z=trans_prob.tolist(),
x=STATE_NAMES,
y=STATE_NAMES,
colorscale="Blues",
zmin=0,
zmax=1,
text=[[f"{v:.0%}" for v in row] for row in trans_prob],
texttemplate="%{text}",
textfont=dict(size=20),
hovertemplate="today: %{y}<br>tomorrow: %{x}<br>probability: %{z:.1%}<extra></extra>",
showscale=False,
)
)
fig.update_layout(
title="Real transition matrix, fit on 4 years of Seattle weather",
xaxis_title="Tomorrow",
yaxis_title="Today",
)
fig.update_yaxes(autorange="reversed")
show(fig, "transition-heatmap", height=380)
Reading the matrix¶
A dry day is followed by another dry day about 76% of the time; a wet day is followed by another wet day about 67% of the time. Both states are "sticky" — weather tends to persist — but dry spells are slightly more self-reinforcing than wet ones, at least in Seattle.
Here's the same information as a flow: out of every 1,000 days in each state, where does the chain send them tomorrow?
fig = go.Figure(
data=go.Sankey(
node=dict(
label=["Dry (today)", "Wet (today)", "Dry (tomorrow)", "Wet (tomorrow)"],
color=["#c98a2b", "#1b6ca8", "#c98a2b", "#1b6ca8"],
pad=25,
thickness=18,
),
link=dict(
source=[0, 0, 1, 1],
target=[2, 3, 2, 3],
value=(trans_prob.flatten() * 1000).tolist(),
color=[
"rgba(201,138,43,0.4)",
"rgba(201,138,43,0.4)",
"rgba(27,108,168,0.4)",
"rgba(27,108,168,0.4)",
],
hovertemplate="%{value:.0f} of every 1000 days<extra></extra>",
),
)
)
fig.update_layout(
title="Where does today's weather send tomorrow's? (per 1,000 days)", font_size=12
)
show(fig, "transition-sankey", height=380)
Does the chain even agree with itself?¶
If we let this chain run forever, the fraction of days it spends in each state should settle into a fixed stationary distribution — and that should match the actual fraction of wet/dry days we observed. If it didn't, something would be wrong with the model. Let's check.
eigvals, eigvecs = np.linalg.eig(trans_prob.T)
stationary = eigvecs[:, np.isclose(eigvals, 1)].flatten().real
stationary = stationary / stationary.sum()
observed = np.array([1 - states.mean(), states.mean()])
comparison = pd.DataFrame(
{
"Chain's stationary distribution": stationary,
"Actually observed in the data": observed,
},
index=STATE_NAMES,
).round(3)
comparison
| Chain's stationary distribution | Actually observed in the data | |
|---|---|---|
| Dry | 0.573 | 0.574 |
| Wet | 0.427 | 0.426 |
They match almost exactly — the chain's long-run behavior is internally consistent with the real climate of Seattle. That's reassuring, but it's a weak test: any correctly-fit Markov chain will pass it by construction. The real test is whether the memoryless assumption itself holds.
The real test: does yesterday actually matter?¶
The Markov property says $P(\text{tomorrow} \mid \text{today}) = P(\text{tomorrow} \mid \text{today}, \text{yesterday}, \ldots)$ — knowing more history shouldn't change the forecast. Let's check that directly: split the data by both today's and yesterday's state, and see whether the forecast for tomorrow actually shifts.
yesterday, today2, tomorrow2 = states[:-2], states[1:-1], states[2:]
rows = []
for today_state, today_label in enumerate(STATE_NAMES):
p_markov = tomorrow2[today2 == today_state].mean()
for yday_state, yday_label in enumerate(STATE_NAMES):
mask = (today2 == today_state) & (yesterday == yday_state)
p_full_history = tomorrow2[mask].mean()
rows.append(
{
"today": today_label,
"yesterday": yday_label,
"n_days": int(mask.sum()),
"P(wet tomorrow | today only)": p_markov,
"P(wet tomorrow | today AND yesterday)": p_full_history,
}
)
memory_test = pd.DataFrame(rows)
memory_test.round(3)
| today | yesterday | n_days | P(wet tomorrow | today only) | P(wet tomorrow | today AND yesterday) | |
|---|---|---|---|---|---|
| 0 | Dry | Dry | 632 | 0.243 | 0.187 |
| 1 | Dry | Wet | 204 | 0.243 | 0.417 |
| 2 | Wet | Dry | 204 | 0.673 | 0.657 |
| 3 | Wet | Wet | 419 | 0.673 | 0.680 |
fig = go.Figure()
for today_state, today_label in enumerate(STATE_NAMES):
subset = memory_test[memory_test["today"] == today_label]
fig.add_trace(
go.Bar(
name=f"today = {today_label}",
x=[f"yesterday = {y}" for y in subset["yesterday"]],
y=(subset["P(wet tomorrow | today AND yesterday)"] * 100).tolist(),
hovertemplate="%{x}<br>P(wet tomorrow)=%{y:.1f}%<extra></extra>",
)
)
fig.add_hline(
y=subset["P(wet tomorrow | today only)"].iloc[0] * 100,
line=dict(dash="dash", color="crimson", width=1.5),
annotation_text=f"Markov forecast for today={today_label}: {subset['P(wet tomorrow | today only)'].iloc[0]:.0%}",
annotation_position="top left" if today_state == 0 else "bottom left",
)
fig.update_layout(
title="If yesterday truly didn't matter, every bar would sit on its dashed line",
yaxis_title="P(wet tomorrow), %",
barmode="group",
)
show(fig, "memory-test", height=460)
So... does it matter?¶
Mostly not, with one honest exception. When today is wet, tomorrow's forecast barely moves whether yesterday was wet (68%) or dry (66%) — the Markov assumption holds up well here. But when today is dry, yesterday matters a lot: if yesterday was also dry, an 81% chance of staying dry tomorrow; if yesterday was wet, that drops to 58%. A dry day sandwiched right after rain is meaningfully less "settled" than a dry day in the middle of a dry spell.
So the memoryless assumption isn't exactly true — there's a real, measurable echo of yesterday, especially around the tail end of a wet spell. But it's a second-order effect, not a first-order one. For a first pass at modeling day-to-day persistence, "today" is doing most of the work.
What the chain gets right: the long-run statistics¶
Even an imperfect chain can be very useful if it reproduces the right long-run behavior. Let's simulate 40 independent weather histories from our fitted 2-state chain and watch the running fraction of wet days settle down — the same "many plausible futures" idea from the ensemble forecasts in the last post, just with a much simpler engine underneath.
rng = np.random.default_rng(11)
n_chains, n_days = 40, 150
start_state = 0 # start dry
running_wet_frac = np.zeros((n_chains, n_days))
for c in range(n_chains):
state = start_state
wet_count = 0
for d in range(n_days):
state = rng.choice(2, p=trans_prob[state])
wet_count += state
running_wet_frac[c, d] = wet_count / (d + 1)
fig = go.Figure()
days = list(range(1, n_days + 1))
for c in range(n_chains):
fig.add_trace(
go.Scatter(
x=days,
y=running_wet_frac[c].tolist(),
mode="lines",
line=dict(color="#1b6ca8", width=1),
opacity=0.25,
hoverinfo="skip",
showlegend=False,
)
)
fig.add_trace(
go.Scatter(
x=days,
y=running_wet_frac.mean(axis=0).tolist(),
mode="lines",
line=dict(color="crimson", width=2.5),
name="Average across simulations",
)
)
fig.add_hline(
y=float(observed[1]),
line=dict(dash="dash", color="black", width=1.5),
annotation_text=f"Actually observed in Seattle: {observed[1]:.1%} of days",
)
fig.update_layout(
title="40 simulated weather histories, all starting dry",
xaxis_title="Day of simulation",
yaxis_title="Running fraction of wet days",
legend=dict(orientation="h", y=1.12),
)
show(fig, "mc-convergence", height=460)
Every simulated run wobbles wildly for the first couple of weeks — pure luck in a short run — but they all settle in around the same long-run wet fraction that we actually observed in the real data. This is the payoff of the Markov chain: even a slightly-wrong, memoryless model can nail the climatology (the long-run statistics) while being genuinely uncertain about any individual day. That's the exact same lesson as the ensemble forecasts from last time, at a fraction of the computational cost.
Why "good enough" works here¶
Real weather systems — fronts, high- and low-pressure systems, atmospheric rivers — tend to move and evolve over a few days, not a few hours. That means today's weather already encodes most of the recent history implicitly: if it's raining today, there's a good chance a system has been sitting overhead for a day or two already, whether or not we explicitly told the model that. A first-order Markov chain isn't ignoring memory so much as absorbing it into "today's state."
In terms of the last post: a Markov chain doesn't fight the atmosphere's chaos at all. It just harvests the short (1–3 day) stretch of genuine persistence that survives before chaos erases it, using the cheapest possible model that can do that.
Where it breaks down¶
Our memory test already caught the crack: dry days that follow a wet spell behave differently from dry days deep in a dry spell. A first-order chain misses that. Two common fixes:
- Higher-order Markov chains — condition on the last k days instead of just 1, trading simplicity for a bit more accuracy (and needing a lot more data to estimate reliably).
- Slower, non-Markov memory — the real long-range predictability in weather doesn't come from yesterday at all, but from boundary conditions that persist for months: ocean temperatures, El Niño/La Niña, soil moisture. That's exactly the seasonal-outlook machinery from the last post, and it lives on a completely different timescale than anything a short-lag Markov chain can see.
The takeaway¶
- A Markov chain assumes tomorrow depends only on today — a bold, "memoryless" simplification.
- Fit on four real years of Seattle weather, the resulting 2-state chain is internally consistent: its long-run stationary distribution matches the actual observed climate.
- The memoryless assumption is mostly true but not perfectly true — yesterday leaves a real, measurable fingerprint, especially on how a wet spell tapers off into dry days.
- Even so, simulated chains reproduce the right long-run statistics, the same payoff we saw from full ensemble weather forecasts, at a tiny fraction of the cost.
- It works because day-to-day weather systems already persist for a few days on their own — the chain isn't defeating chaos, it's cheaply capturing the short window of persistence that chaos hasn't erased yet.