twin ridges: visualizing housing expensiveness trends across US states

Housing affordability is a constant topic of conversation in the United States, largely because affordability is so hard to find. The ratio of home values to household income provides a simple way to measure how expensive it is to buy a home relative to what typical households can afford. But when you look at this ratio across all 50 states over the past 25 years, a surprising pattern emerges: despite their differences, nearly every state follows a similar “twin-peaked” pattern of housing expensiveness, with peaks occurring around two distinct periods— the pre-2008 housing bubble and the post-COVID housing market surge.

In this analysis, I’ll combine Zillow’s Home Value Index (ZHVI) with median household income data from FRED to create a ridgeline plot that visualizes these trends across all states. The visualization reveals that while 26 states reached their peak expensiveness in 2022, 15 states actually peaked before the 2008 recession, suggesting that income growth in some regions— particularly in the Northeast— has partially offset home price appreciation in recent years relative to the pre-2008 peak.

This notebook walks through the data collection, processing, and visualization steps needed to create this analysis, including some methodological choices that warrant discussion along the way.

This post is another attempt at visualizing the same data that I covered in my previous post on housing expensiveness; the visualization there tried to cram too much information into a single plot, whereas this version focuses on the twin-peaked pattern across all states.

libraries and imports

This is where we’ll import the python packages used for this analysis. Reach out if you need installation details for package versions, or even a copy of the notebook file.

from pathlib import Path
import pandas as pd
import polars as pl
# data fetching
import pandas_datareader.data as web # FRED downloads
import requests
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import dates as mdates
from matplotlib.lines import Line2D
# transforms are used for text positioning
from matplotlib import transforms

# customize font in visualizations
mpl.rcParams['font.family'] = 'monospace'
mpl.rcParams['font.size'] = 14
mpl.rcParams['figure.facecolor'] = '#fafaf9'

labeling states with abbreviations

The Zillow data we’ll download shortly uses state names, not USPS codes; FRED uses USPS codes in its series tickers, hence defining them before that download. We’ll also use these 2-letter codes to label subplots in the grid cartogram.

state_to_abbr = { 'Alabama': 'AL', 'Alaska': 'AK', 'Arizona': 'AZ', 'Arkansas': 'AR', 'California': 'CA', 'Colorado': 'CO', 'Connecticut': 'CT', 'Delaware': 'DE', 'Florida': 'FL', 'Georgia': 'GA', 'Hawaii': 'HI', 'Idaho': 'ID', 'Illinois': 'IL', 'Indiana': 'IN', 'Iowa': 'IA', 'Kansas': 'KS', 'Kentucky': 'KY', 'Louisiana': 'LA', 'Maine': 'ME', 'Maryland': 'MD', 'Massachusetts': 'MA', 'Michigan': 'MI', 'Minnesota': 'MN', 'Mississippi': 'MS', 'Missouri': 'MO', 'Montana': 'MT', 'Nebraska': 'NE', 'Nevada': 'NV', 'New Hampshire': 'NH', 'New Jersey': 'NJ', 'New Mexico': 'NM', 'New York': 'NY', 'North Carolina': 'NC', 'North Dakota': 'ND', 'Ohio': 'OH', 'Oklahoma': 'OK', 'Oregon': 'OR', 'Pennsylvania': 'PA', 'Rhode Island': 'RI', 'South Carolina': 'SC', 'South Dakota': 'SD', 'Tennessee': 'TN', 'Texas': 'TX', 'Utah': 'UT', 'Vermont': 'VT', 'Virginia': 'VA', 'Washington': 'WA', 'West Virginia': 'WV', 'Wisconsin': 'WI', 'Wyoming': 'WY', 'District of Columbia': 'DC', 'United States': ''} # we'll fill USA in later, this needs to be empty for FRED ticker to match

state_to_abbr = pd.DataFrame.from_dict(state_to_abbr, orient='index')
state_to_abbr.columns = ['usps']
state_to_abbr.index.name = 'state'
state_to_abbr.reset_index(inplace=True)

home value data: Zillow Home Value Index (ZHVI)

Both Zillow and Redfin publish monthly home value indices, but Zillow’s goes farther back— 2000 vs 2012— so that’s what I’ll be using. Zillow’s Home Value Index (ZHVI) is a smoothed estimate of the median home value in a given geography, calculated by taking the average of the middle third of the estimated home values in a given month.

The ZHVI is based on estimates, rather than actual transactions, which might draw skepticism. This is a concern that Zillow’s economic research team has gone to lengths to address: Zillow’s research articles explaining the ZHVI point out that the estimation errors for houses aren’t systematically biased when compared to later realized sale prices, implying that an index built from their estimated home values aggregated over a sufficiently large number of properties should be an unbiased estimator of home value.

And realistically, any index that claims to measure the “median home value in a state” is going to rely on a slew of simplifying assumptions. This data passes a smell test, at least to my nose.

def download_if_not_exists(url, save_path: Path | None = None):
if save_path is None:
# assume saving to current directory
save_path = Path(url).name
if not Path(save_path).exists():
response = requests.get(url)
with open(save_path, 'wb') as f:
f.write(response.content)
return Path(save_path).absolute()

# note: to fetch updated data, delete the CSVs to trigger re-download.
# ZHVI updates monthly on the 16th. Automating data refresh checking is out of scope here.
zhvi_national_path = download_if_not_exists("https://files.zillowstatic.com/research/public_csvs/zhvi/Metro_zhvi_uc_sfrcondo_tier_0.33_0.67_sm_sa_month.csv", Path('zhvi_national_data.csv'))
zhvi_national_df = pd.read_csv(zhvi_national_path)

zhvi_path = download_if_not_exists("https://files.zillowstatic.com/research/public_csvs/zhvi/State_zhvi_uc_sfrcondo_tier_0.33_0.67_sm_sa_month.csv", Path('zhvi_data.csv'))
zhvi_df = pd.read_csv(zhvi_path)

# process the ZHVI data
def process_zhvi_data(zhvi_df: pd.DataFrame) -> pd.DataFrame:
zhvi_df = zhvi_df\
.drop(columns=['RegionID', 'SizeRank', 'RegionType', 'StateName'])\
.melt(id_vars=['RegionName'])\
.rename(columns={'RegionName': 'state', 'value': 'zhvi', 'variable': 'date'})\
.assign(date=lambda x: pd.to_datetime(x['date']))\
.set_index(['state', 'date'])\
.sort_index()
return zhvi_df

# filter to national before processing because same file includes metro areas
zhvi_national_df = process_zhvi_data(zhvi_national_df.loc[zhvi_national_df['RegionName'] == 'United States'])

# process state-level ZHVI data
zhvi_df = process_zhvi_data(zhvi_df)

# concatenate national and state-level ZHVI data
zhvi_df = pd.concat([zhvi_national_df, zhvi_df])

income data: median household income from FRED

FRED, the Federal Reserve Bank of St. Louis’s data repository, publishes data from the US Census Bureau on Median Household Income by state. It’s annual data, and the series notes state that household data are collected as of the end of March in each year. Since ZHVI is a monthly series, we can linearly interpolate the annual household income data to provide a more detailed time series.

# template for median household income by state, current dollars, not seasonally adjusted
# includes national data labeled as 'United States'
state_income_tickers = [
f'MEHOINUS{abbr}A646N' for abbr in state_to_abbr['usps'].values
]

state_to_abbr['fred_ticker'] = state_income_tickers

mhhi_path = Path('mhhi_data.csv')
if not mhhi_path.exists():
# annual data are indexed with the first day of the year
hhi_data = web.DataReader(state_income_tickers, 'fred', start='1999-01-01')
hhi_data.to_csv(mhhi_path)
else:
hhi_data = pd.read_csv(mhhi_path, index_col=0, parse_dates=True)

# Step 1: Change dates from 01-01-YYYY to 03-31-YYYY
# Replace the year-start dates with March 31st of each year
hhi_data.index = hhi_data.index.map(
lambda x: pd.Timestamp(year=x.year, month=3, day=31)
)
# Step 2: Convert to monthly data with linear interpolation
# Resample to monthly frequency ('MS' = month start or 'M' = month end) (ME on newer versions of pandas)
# We'll use 'M' for month-end since dates are end-of-month (03-31)
hhi_data = hhi_data.resample('ME').asfreq()
# Step 3: Linearly interpolate all columns
hhi_data = hhi_data.interpolate(method='linear')

hhi_data = hhi_data\
.melt(ignore_index=False)\
.reset_index()\
.rename(columns={'DATE': 'date'})\
.merge(state_to_abbr, left_on='variable', right_on='fred_ticker', how='left')\
.set_index(['state', 'date'])\
.rename(columns={'value': 'median_household_income'})\
.drop(columns=['fred_ticker', 'variable'])

# replace empty USPS with 'USA' (for labeling later)
hhi_data.loc[hhi_data['usps'] == '', 'usps'] = 'USA'

combining the home value and income data

Here’s a modeling assumption that I’m sure some are going to find controversial: I’m going to extrapolate the latest valid data point (corresponding to 3/31/2024) forward using each state’s average growth rate over the period 2000-2024.

All states show a general upward trend in median household income, if only because of inflation, and growth at the average rate achieved from 2000-2024 seems like a reasonable baseline assumption, even if it does feel like we’re living through extraordinary economic times. At the very least, this is a more reasonable assumption than forward-filling the income data, since it’s likely that household income has at least partially kept up with inflation.

state_data = zhvi_df.merge(hhi_data, left_index=True, right_index=True, how='left').sort_index()
# forward-fill the USPS code (not controversial)
state_data['usps'] = state_data.groupby(level=0)['usps'].ffill()
# first linearly interpolate missing values in the ZHVI data
# (some periods have missing data due to insufficient transaction volume/data collection issues)
state_data['zhvi'] = state_data.groupby(level=0)['zhvi'].apply(lambda x: x.interpolate(method='linear', limit_direction='forward', limit_area='inside')).droplevel(0, 'index')

def extrapolate_series(series: pd.Series) -> pd.Series:
"""
Calculate average growth rate from non-null values and extrapolate forward
from the last known value to avoid discontinuities.
Verbose boilerplate, by Claude.
"""
# Get the non-null values
non_null_mask = series.notna()

if non_null_mask.sum() < 2:
# Not enough data points to calculate growth rate
return series

# Get non-null values
non_null_values = series[non_null_mask]

# Calculate period-over-period growth rates
growth_rates = non_null_values.pct_change().dropna()

if len(growth_rates) == 0:
# Can't calculate growth rate
return series

# Calculate average growth rate
avg_growth_rate = growth_rates.mean()

# Get the last known value and its position
last_known_idx = non_null_values.index[-1]
last_known_value = non_null_values.iloc[-1]

# Find positions after the last known value
last_position = series.index.get_loc(last_known_idx)

# Fill forward with compound growth
result = series.copy()
current_value = last_known_value

for i in range(last_position + 1, len(series)):
if pd.isna(result.iloc[i]):
current_value = current_value * (1 + avg_growth_rate)
result.iloc[i] = current_value

return result

# Apply the extrapolation
state_data['median_household_income'] = (
state_data.groupby(level=0, group_keys=False)['median_household_income']
.apply(extrapolate_series)
)

state_data['expensiveness_index'] = state_data['zhvi'] / state_data['median_household_income']

ridges and trends

Canonically, “ridgeline plots” like the one I’m making in this post are used to show how the distribution of a variable changes over time or across groups. Because they don’t include a scale for the Y-axis, you shouldn’t use a ridgeline plot to show a time series that you want the reader to be able to pick out specific values from.

In this case, however, the point I’m trying to make is less about the exact “housing expensiveness” values each state reached than it’s about broad patterns in the timing of peaks and troughs across states. I was surprised to find that, despite their differences, nearly every state followed a similar “twin-peaked” housing expensiveness pattern over the past 25 years.

I also found it interesting that housing expensiveness in many states peaked before the 2008 recession. Despite significant home price appreciation in almost every state since the COVID-19 pandemic, many expensive housing markets, particularly in the Northeast, have seen median household incomes grow sufficiently that the ratio of home value to household income remains below pre-2008 peaks.

”excess expensiveness”

Some states have housing markets that are persistently more expensive than other parts of the country— for example, Hawaii, California, and Massachusetts. There’s a certain baseline expensiveness that’s structurally built into each state’s housing market, driven by underlying factors like the ease of building new housing. In order to highlight common trends across states, I’m going to compare observed housing expensiveness ratio values to the minimum observed ratio for each state. In any other treatment of the data, differences in minimum values across states conceal shared responses to national trends.

Here’s the data manipulation code for the “excess expensiveness” metric:

# computing "excess expensiveness" as a percent (decimal, 0-1) of each state's minimum observed ratio
# switch to polars; not mapping, don't need pandas.
sd = pl.from_pandas(state_data.reset_index())
sd = sd.with_columns(
pl.col('expensiveness_index').min().over('state').alias('min_expensiveness')
).with_columns(
(pl.col('expensiveness_index') - pl.col('min_expensiveness')).alias('excess_expensiveness')
).with_columns(
(pl.col('excess_expensiveness') / pl.col('min_expensiveness')).alias('excess_expensiveness_pc'),
pl.col('excess_expensiveness').max().over('state').alias('max_excess_expensiveness')
).with_columns(
pl.col('excess_expensiveness_pc').max().over('state').alias('max_excess_expensiveness_pc'),
pl.col('date').filter(pl.col('excess_expensiveness') == pl.col('max_excess_expensiveness')).first().over('state').alias('date_of_max_excess_expensiveness')
).sort(['state', 'date'], descending=[False, False])

number of states with expensiveness peaking in a given year

Total adds up to 51 because this data includes the District of Columbia.

sd.filter(pl.col('state') != 'United States').with_columns(
pl.col('date_of_max_excess_expensiveness').dt.year().alias('year_of_max_excess_expensiveness')
).group_by('year_of_max_excess_expensiveness')\
.agg(pl.col('state').n_unique().alias('n_states'))\
.sort('year_of_max_excess_expensiveness', descending=True)\
.style\
.cols_label({
'year_of_max_excess_expensiveness': 'Year of Peak Expensiveness',
'n_states': 'Number of States'
})\
.data_color('n_states', palette='Blues')\
.tab_source_note('Source: Zillow Home Value Index (zillow.com/research/data), FRED (fred.stlouisfed.org), aaronjbecker.com')\
.tab_header(f'Number of States with Expensiveness Peaking in a Given Year',
subtitle='Expensiveness measured as Median Home Value / Median Household Income')

Year of Peak Expensiveness	Number of States
Number of States with Expensiveness Peaking in a Given Year
Expensiveness measured as Median Home Value / Median Household Income
2025	5
2024	1
2022	26
2021	1
2016	1
2015	1
2010	1
2008	1
2007	2
2006	7
2005	4
2004	1
Source: Zillow Home Value Index (zillow.com/research/data), FRED (fred.stlouisfed.org), aaronjbecker.com

creating the ridgeline plot

If this matplotlib plotting code seems verbose and complicated to you, you’re right! Building almost any non-standard visualization with matplotlib requires diving into its thorny imperative layout processes. Fortunately, AI coding assistants like Claude and Cursor can do a decent job generating code to your specifications, although you’ll often have to nudge them in the right direction.

I still glance at R’s ggplot2 community with envy from time to time, but then I remember that every abstraction ends up leaking. To some extent LLMs have become my matplotlib abstraction layer; since I still understand the underlying code, I can jump in and tweak things to my heart’s content.

# Declare variables for the series and color mapping
series_variable = 'excess_expensiveness_pc'
# used for peak value marking
peak_variable = 'max_excess_expensiveness_pc'
label_variable = 'expensiveness_index'

even_color = '#155dfc'
odd_color = '#51a2ff'
usa_color = '#fb2c36'
grid_color = '#79716b'
# color of vertical line and filled circle at max value
max_color = '#ffd230'
# background color for figure
face_color = '#e7e5e4'

# Convert to pandas for easier matplotlib integration
df = sd.to_pandas()

# Get unique states, sorted by metric selected above
states = sd.select(pl.col('state').unique(maintain_order=True)).to_series().to_list()

# Number of subplots
n_states = len(states)
height = 1.0

# Create figure with subplots
fig, axes = plt.subplots(n_states, 1, figsize=(15, height * n_states),
facecolor=(0, 0, 0, 0), sharex=True, sharey=False)

fig.set_facecolor(face_color)

# positions for the min/max/last labels (in subplot axis coordinates)
min_pos_x = 0.975
max_pos_x = 1.035
last_pos_x = 1.095

# Plot each state
for idx, state in enumerate(states):
ax = axes[idx]

# Get data for this state
state_data = df[df['state'] == state].sort_values('date')
usps = state_data['usps'].iloc[0]
dates = state_data['date'].values
expensiveness = state_data[series_variable].values
peak_val = state_data[peak_variable].iloc[0]

# Find the maximum value and corresponding date for this state
max_idx = state_data[series_variable].idxmax()
max_value = state_data.loc[max_idx, series_variable]
max_date = state_data.loc[max_idx, 'date']

# Find the minimum value and corresponding date for this state
min_idx = state_data[series_variable].idxmin()
min_date = state_data.loc[min_idx, 'date']

# Get min and max values of label_variable for this state
min_label_value = state_data[label_variable].min()
max_label_value = state_data[label_variable].max()
last_label_value = state_data[label_variable].iloc[-1]

# Get color for this state: USA gets usa_color, others alternate between even/odd
if usps == 'USA':
color = usa_color
elif idx % 2 == 0:
color = even_color
else:
color = odd_color

# Fill the area under the curve with the color
# Set zorder to ensure it appears above figure-level grid (zorder=-10)
ax.fill_between(dates, 0, expensiveness, color=color, alpha=0.7, zorder=1)

# Plot the line in black
ax.plot(dates, expensiveness, color='black', linewidth=2, zorder=2)

# Plot a thin black vertical line from axis to max value
ax.plot([max_date, max_date], [0, max_value], color=max_color, linewidth=2.5, zorder=3)

# Plot a heavy black mark on the x-axis at the minimum value date
ax.scatter([min_date], [0], marker='o', s=150, color='black', facecolor='none', zorder=5, clip_on=False)

ax.scatter([max_date], [max_value], marker='o', s=150, color=max_color, zorder=5, clip_on=False)

# Add horizontal reference line at y=0, extended to run under labels including Range column
# Use blended transform: x in axes coordinates, y in data coordinates
blended_transform = transforms.blended_transform_factory(ax.transAxes, ax.transData)
ax.plot([0, 1.15], [0, 0], linewidth=2, linestyle='-', color=color,
clip_on=False, transform=blended_transform)

# Add state label to the left
ax.text(-0.01, 0.2, usps, fontweight='bold', color=color, fontsize=30,
ha='right', va='center', transform=ax.transAxes)

# Add table-like labels to the right
# Column headers only on the first subplot
value_label_font_size = 16
if idx == 0:
ax.text(min_pos_x, 1.05, 'Min', fontweight='bold', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=20)
ax.text(max_pos_x, 1.05, 'Max', fontweight='bold', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=20)
ax.text(last_pos_x, 1.05, 'Last', fontweight='bold', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=20)
# Values for all subplots, moved to the left
ax.text(min_pos_x, 0.2, f'{min_label_value:.1f}x', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=value_label_font_size)
ax.text(max_pos_x, 0.2, f'{max_label_value:.1f}x', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=value_label_font_size)
ax.text(last_pos_x, 0.2, f'{last_label_value:.1f}x', color='black',
ha='left', va='center', transform=ax.transAxes, fontsize=value_label_font_size)

# Remove unnecessary axes details
ax.set_yticks([])
ax.set_ylabel('')
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

# Set transparent background
ax.set_facecolor((0, 0, 0, 0))

# Hide x-axis tick marks (grid will be drawn on figure level)
ax.tick_params(axis='x', which='major', length=0, width=0, labelsize=0)

# Set the subplots to overlap
fig.subplots_adjust(hspace=-0.5)

###############
# GRID LABELS #
###############
# Get x-axis tick locations from any axis (they're all shared)
reference_ax = axes[-1]  # Use last axis to ensure ticks are set
xticks = reference_ax.get_xticks()
# Filter to valid tick locations within the view
xlim = reference_ax.get_xlim()
xtick_locations = [tick for tick in xticks if pd.notna(tick) and xlim[0] <= tick <= xlim[1]]

# Convert matplotlib numeric tick locations to dates using matplotlib's date utilities
xtick_dates = [mdates.num2date(loc) for loc in xtick_locations]
xtick_years = [str(dt.year) for dt in xtick_dates]

# Get bounding boxes for positioning labels
top_ax = axes[0]
bottom_ax = axes[-1]

# Get the position of the top and bottom axes in figure coordinates
top_bbox = top_ax.get_position()  # Returns Bbox with x0, y0, x1, y1 in figure coordinates
bottom_bbox = bottom_ax.get_position()

# Position labels just above the top subplot and below the bottom subplot
# Add a small offset (0.02 in figure coordinates) above/below the subplots
top_y = top_bbox.y1  # Just above the top subplot
bottom_y = bottom_bbox.y0  # Just below the bottom subplot

# Create blended transforms for positioning labels
# x in data coordinates, y in figure coordinates
top_transform = transforms.blended_transform_factory(
top_ax.transData,  # x in data coordinates from first axis
fig.transFigure     # y in figure coordinates (0 to 1)
)
bottom_transform = transforms.blended_transform_factory(
bottom_ax.transData,  # x in data coordinates from last axis
fig.transFigure       # y in figure coordinates (0 to 1)
)

# Add year labels at the top (above first subplot)
top_label_artists = []
for tick_loc, year in zip(xtick_locations, xtick_years):
# Position just above the top subplot
text_artist = fig.text(tick_loc, top_y, year, ha='center', va='bottom',
transform=top_transform, fontsize=20, color='black')
top_label_artists.append(text_artist)

# Add year labels at the bottom (below last subplot)
bottom_label_artists = []
for tick_loc, year in zip(xtick_locations, xtick_years):
# Position just below the bottom subplot
text_artist = fig.text(tick_loc, bottom_y, year, ha='center', va='top',
transform=bottom_transform, fontsize=20, color='black')
bottom_label_artists.append(text_artist)

##############
# GRID LINES #
##############
# Draw grid lines on the figure level, from bottom of top labels to top of bottom labels
# Get the actual bbox of a label to determine label height
# Need to render the figure first to get accurate bboxes
fig.canvas.draw()
# we can use these locations to position text as well
grid_top_y = None
grid_bottom_y = None
if top_label_artists:
# Get bbox of top label in display coordinates (pixels)
renderer = fig.canvas.get_renderer()
top_label_bbox_display = top_label_artists[0].get_window_extent(renderer=renderer)
# Convert bbox to figure coordinates using the figure's transform
# The bbox is in display coordinates, we need figure coordinates
fig_bbox_display = fig.get_window_extent()
# Display coords: y0 is bottom, y1 is top (both in pixels from figure bottom)
# Figure coords: 0 is bottom, 1 is top (normalized)
top_label_bottom_display = top_label_bbox_display.y0
top_label_bottom_fig = (top_label_bottom_display - fig_bbox_display.y0) / fig_bbox_display.height
grid_top_y = top_label_bottom_fig  # Just below the bottom of top labels
else:
# Fallback if no labels
grid_top_y = top_y - 0.015

if bottom_label_artists:
# Get bbox of bottom label in display coordinates (pixels)
renderer = fig.canvas.get_renderer()
bottom_label_bbox_display = bottom_label_artists[0].get_window_extent(renderer=renderer)
# Convert bbox to figure coordinates
fig_bbox_display = fig.get_window_extent()
# Get the top of the bottom label bbox in figure coordinates
bottom_label_top_display = bottom_label_bbox_display.y1
bottom_label_top_fig = (bottom_label_top_display - fig_bbox_display.y0) / fig_bbox_display.height
grid_bottom_y = bottom_label_top_fig  # Just above the top of bottom labels
else:
# Fallback if no labels
grid_bottom_y = bottom_y + 0.015

# Create a blended transform for drawing grid lines: x in data coordinates, y in figure coordinates
grid_transform = transforms.blended_transform_factory(
reference_ax.transData,  # x in data coordinates
fig.transFigure           # y in figure coordinates
)

# Draw vertical grid lines at each tick location using Line2D objects
for tick_loc in xtick_locations:
line = Line2D([tick_loc, tick_loc], [grid_bottom_y, grid_top_y],
color=grid_color, linewidth=1, alpha=1, zorder=-10,
transform=grid_transform, clip_on=False)
fig.add_artist(line)

####################
# TABLE GRID LINES #
####################
# draw lines between min/max/last
table_col_transform = transforms.blended_transform_factory(
reference_ax.transAxes,
fig.transFigure
)

# draw line between min/max
col_position_offset = 0.0075
line = Line2D([max_pos_x - col_position_offset, max_pos_x - col_position_offset], [grid_bottom_y, grid_top_y], color='black', linewidth=1, transform=table_col_transform)
fig.add_artist(line)

# draw line between max/last
line = Line2D([last_pos_x - col_position_offset, last_pos_x - col_position_offset], [grid_bottom_y, grid_top_y], color='black', linewidth=1, transform=table_col_transform)
fig.add_artist(line)

#######################
# TOP AND BOTTOM TEXT #
#######################
mainTitle = 'Timing of Peak Housing Expensiveness\n by US State, Jan. 2000 - Sept. 2025'
subtitle = 'Change in (Median Home Price / Median Household Income) vs. State Min.'
subtitle2 = 'States are alphabetically. Each y-axis is scaled to that state\'s range.'
subtitle3 = 'Vertical lines indicate date of max expensiveness. Empty circles indicate date of minimum expensiveness.'
subtitle4 = '26 states were most expensive in 2022, but 15 peaked prior to the 2008 recession.'

# Top text
text = fig.text(0.05, grid_top_y + 0.01, subtitle4, fontsize=18, ha='left', fontfamily='sans-serif', color='#44403b', fontweight='bold', va='bottom')
ex = text.get_window_extent()
x, y = text.get_position()
t = transforms.offset_copy(text._transform, y=ex.height + 10, units='dots')
text = fig.text(x, y, subtitle3, transform=t, va='bottom', fontfamily='sans-serif', fontsize=18, color='#44403b')
ex = text.get_window_extent()
x, y = text.get_position()
t = transforms.offset_copy(text._transform, y=ex.height + 10, units='dots')
text = fig.text(x, y, subtitle2, transform=t, va='bottom', fontfamily='sans-serif', fontsize=18, color='#44403b')
ex = text.get_window_extent()
x, y = text.get_position()
t = transforms.offset_copy(text._transform, y=ex.height + 20, units='dots')
text = fig.text(x, y, subtitle, transform=t, va='bottom', fontsize=24)
ex = text.get_window_extent()
x, y = text.get_position()
t = transforms.offset_copy(text._transform, y=ex.height + 10, units='dots')
text = fig.text(x, y, mainTitle, transform=t, va='bottom', fontsize=44, fontweight='bold')

# Bottom text
text = fig.text(0.05, grid_bottom_y - 0.0125,
"Note: ND, MT, NM, and WY are missing early ZHVI data. Median Household Income from 3/31/2024 - 9/30/2025\n   is extrapolated using each state's 25Y trend growth rate.", fontsize=18, ha='left', fontfamily='sans-serif', color='#44403b', va='top')
ex = text.get_window_extent()
x, y = text.get_position()
t = transforms.offset_copy(text._transform, y=-ex.height - 10, units='dots')
text = fig.text(x, y, "Source: Zillow Home Value Index (zillow.com/research/data),\n    FRED (fred.stlouisfed.org)", transform=t, va='top', fontfamily='sans-serif', fontsize=18, color='#44403b')

# website link (we can re-use the y position and transform from source note)
text = fig.text(1.01, y,
"aaronjbecker.com", fontsize=36, ha='right', va='top', transform=t)

examining the trends

Why did so many states, particularly in the Northeast, experience housing expensiveness peaks in the pre-2008 housing market cycle instead of the current post-COVID cycle? We can look to the numerator and denominator of the ratio we’re using to measure expensiveness for possible explanations: either home values increased relatively less than in other states, or incomes increased relatively more.

# Faceted line chart: zhvi and median income by cycle group
from datetime import datetime

# filter states by cycle of peak expensiveness (dates are rough and manually chosen)
covid_peak = ((pl.col('date_of_max_excess_expensiveness') >= datetime(2020, 1, 1)) & (pl.col('date_of_max_excess_expensiveness') <= datetime(2024, 12, 31)))
pre_2008_peak = ((pl.col('date_of_max_excess_expensiveness') >= datetime(2004, 1, 1)) & (pl.col('date_of_max_excess_expensiveness') <= datetime(2008, 12, 31)))

# Get unique states for each cycle group
covid_states = sd.filter(covid_peak).select('state').unique().to_series().to_list()
pre_2008_states = sd.filter(pre_2008_peak).select('state').unique().to_series().to_list()

# Filter data for each cycle group
sd_covid_full = sd.filter(pl.col('state').is_in(covid_states))
sd_pre_2008_full = sd.filter(pl.col('state').is_in(pre_2008_states))

# Index each state to 100 at its starting value for both zhvi and median household income
# Get the first (earliest) value for each state
sd_covid_indexed = sd_covid_full.sort('state', 'date').with_columns(
pl.col('zhvi').first().over('state').alias('zhvi_start'),
pl.col('median_household_income').first().over('state').alias('income_start')
).with_columns(
(pl.col('zhvi') / pl.col('zhvi_start') * 100).alias('zhvi_indexed'),
(pl.col('median_household_income') / pl.col('income_start') * 100).alias('median_household_income_indexed')
).drop(['zhvi_start', 'income_start'])

sd_pre_2008_indexed = sd_pre_2008_full.sort('state', 'date').with_columns(
pl.col('zhvi').first().over('state').alias('zhvi_start'),
pl.col('median_household_income').first().over('state').alias('income_start')
).with_columns(
(pl.col('zhvi') / pl.col('zhvi_start') * 100).alias('zhvi_indexed'),
(pl.col('median_household_income') / pl.col('income_start') * 100).alias('median_household_income_indexed')
).drop(['zhvi_start', 'income_start'])

# Convert to pandas for easier plotting
df_covid = sd_covid_indexed.to_pandas()
df_pre_2008 = sd_pre_2008_indexed.to_pandas()

# Define colors
covid_color = '#155dfc'
pre_2008_color = '#e7000b'
face_color = '#fafaf9'
grid_color = '#79716b'

# Create 2x2 subplot grid
fig, axes = plt.subplots(2, 2, figsize=(16, 12), facecolor=face_color)
fig.set_facecolor(face_color)

# Top row: zhvi
# Left: pre-2008
ax_zhvi_pre2008 = axes[0, 0]
for state in pre_2008_states:
state_data = df_pre_2008[df_pre_2008['state'] == state].sort_values('date')
ax_zhvi_pre2008.plot(state_data['date'], state_data['zhvi_indexed'],
color=pre_2008_color, alpha=0.4, linewidth=1)

# Calculate mean zhvi for pre-2008 states
mean_zhvi_pre2008 = df_pre_2008.groupby('date')['zhvi_indexed'].mean().reset_index()
ax_zhvi_pre2008.plot(mean_zhvi_pre2008['date'], mean_zhvi_pre2008['zhvi_indexed'],
color=pre_2008_color, linewidth=3, label='Mean', zorder=10)

ax_zhvi_pre2008.set_title('Pre-2008 Peak Cycle', fontsize=16, fontweight='bold', pad=15)
ax_zhvi_pre2008.set_ylabel('ZHVI (Indexed to 100 at Start)', fontsize=14, fontweight='bold')
ax_zhvi_pre2008.grid(True, alpha=0.3, color=grid_color)
ax_zhvi_pre2008.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
# Format y-axis as index values
ax_zhvi_pre2008.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: f'{x:.0f}'))

# Right: COVID
ax_zhvi_covid = axes[0, 1]
for state in covid_states:
state_data = df_covid[df_covid['state'] == state].sort_values('date')
ax_zhvi_covid.plot(state_data['date'], state_data['zhvi_indexed'],
color=covid_color, alpha=0.4, linewidth=1)

# Calculate mean zhvi for COVID states
mean_zhvi_covid = df_covid.groupby('date')['zhvi_indexed'].mean().reset_index()
ax_zhvi_covid.plot(mean_zhvi_covid['date'], mean_zhvi_covid['zhvi_indexed'],
color=covid_color, linewidth=3, label='Mean', zorder=10)

ax_zhvi_covid.set_title('COVID Peak Cycle', fontsize=16, fontweight='bold', pad=15)
ax_zhvi_covid.set_ylabel('ZHVI (Indexed to 100 at Start)', fontsize=14, fontweight='bold')
ax_zhvi_covid.grid(True, alpha=0.3, color=grid_color)
ax_zhvi_covid.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
# Format y-axis as index values
ax_zhvi_covid.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: f'{x:.0f}'))

# Share y-axis for top row - calculate global min/max for zhvi indexed
zhvi_min = min(df_pre_2008['zhvi_indexed'].min(), df_covid['zhvi_indexed'].min())
zhvi_max = max(df_pre_2008['zhvi_indexed'].max(), df_covid['zhvi_indexed'].max())
ax_zhvi_pre2008.set_ylim(zhvi_min, zhvi_max)
ax_zhvi_covid.set_ylim(zhvi_min, zhvi_max)

# Bottom row: median household income
# Left: pre-2008
ax_income_pre2008 = axes[1, 0]
for state in pre_2008_states:
state_data = df_pre_2008[df_pre_2008['state'] == state].sort_values('date')
ax_income_pre2008.plot(state_data['date'], state_data['median_household_income_indexed'],
color=pre_2008_color, alpha=0.4, linewidth=1)

# Calculate mean income for pre-2008 states
mean_income_pre2008 = df_pre_2008.groupby('date')['median_household_income_indexed'].mean().reset_index()
ax_income_pre2008.plot(mean_income_pre2008['date'], mean_income_pre2008['median_household_income_indexed'],
color=pre_2008_color, linewidth=3, label='Mean', zorder=10)

ax_income_pre2008.set_title('Pre-2008 Peak Cycle', fontsize=16, fontweight='bold', pad=15)
ax_income_pre2008.set_ylabel('Median Household Income\n(Indexed to 100 at Start)', fontsize=14, fontweight='bold')
ax_income_pre2008.set_xlabel('Date', fontsize=14, fontweight='bold')
ax_income_pre2008.grid(True, alpha=0.3, color=grid_color)
ax_income_pre2008.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
# Format y-axis as index values
ax_income_pre2008.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: f'{x:.0f}'))

# Right: COVID
ax_income_covid = axes[1, 1]
for state in covid_states:
state_data = df_covid[df_covid['state'] == state].sort_values('date')
ax_income_covid.plot(state_data['date'], state_data['median_household_income_indexed'],
color=covid_color, alpha=0.4, linewidth=1)

# Calculate mean income for COVID states
mean_income_covid = df_covid.groupby('date')['median_household_income_indexed'].mean().reset_index()
ax_income_covid.plot(mean_income_covid['date'], mean_income_covid['median_household_income_indexed'],
color=covid_color, linewidth=3, label='Mean', zorder=10)

ax_income_covid.set_title('COVID Peak Cycle', fontsize=16, fontweight='bold', pad=15)
ax_income_covid.set_ylabel('Median Household Income\n(Indexed to 100 at Start)', fontsize=14, fontweight='bold')
ax_income_covid.set_xlabel('Date', fontsize=14, fontweight='bold')
ax_income_covid.grid(True, alpha=0.3, color=grid_color)
ax_income_covid.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
# Format y-axis as index values
ax_income_covid.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: f'{x:.0f}'))

# Share y-axis for bottom row - calculate global min/max for income indexed
income_min = min(df_pre_2008['median_household_income_indexed'].min(), df_covid['median_household_income_indexed'].min())
income_max = max(df_pre_2008['median_household_income_indexed'].max(), df_covid['median_household_income_indexed'].max())
ax_income_pre2008.set_ylim(income_min, income_max)
ax_income_covid.set_ylim(income_min, income_max)

# Add title and explanatory note
fig.suptitle('ZHVI and Median Household Income by Cycle Group (Indexed to 100 at Start)',
fontsize=20, fontweight='bold', y=0.98)

# Add explanatory note at the top
fig.text(0.5, 0.95, 'Each line represents a state. States are categorized by the cycle in which expensiveness peaked.\nValues are indexed to 100 at each state\'s starting value to show relative growth.',
fontsize=14, ha='center', va='top', fontfamily='sans-serif', color='#44403b', style='italic')

# Add source notes at the bottom left
fig.text(0.01, 0.00, "Source: Zillow Home Value Index (zillow.com/research/data),\n    FRED (fred.stlouisfed.org)",
fontsize=14, ha='left', va='bottom', fontfamily='sans-serif', color='#44403b')

# Add website link
fig.text(0.99, 0.01, "aaronjbecker.com", fontsize=36, ha='right', va='top')

plt.tight_layout(rect=[0, 0.03, 1, 0.97])
plt.show()

The difference in peak expensiveness timing between groups of states that peaked in each cycle seems to be driven more by home value trends than by trends in median household income. When you normalize against starting values, both groups experienced similar growth in median household income. As a group, states whose expensiveness peaked in the pre-2008 cycle experienced sharper increases and steeper declines in home values pre-2008 than states whose expensiveness peaked in the COVID cycle. There are outliers in the COVID cycle group that experienced a sharp boom-and-bust before and after the 2008 financial crisis, but those states later went on to boom even dramatically in the post-COVID cycle.

As to why states where expensiveness peaked earlier showed larger home value increases before 2008… that’s a question that warrants more detailed examination in a future post. No obvious single narrative, like population growth, makes sense given the states in each grouping.

conclusions

The cost of homeownership receives a lot of attention in the US media, and on social media sites like reddit. Owning your own home has long been seen as a key component of the American Dream, one of the milestones that defines adulthood in our culture. Comparing home values to household income is one way of measuring how attainable this goal is for the average American, and there’s some good news in the data: in every state except Wisconsin, home values are off from their peak expensiveness.

The bad news is that the reason for the home value decline, increased mortgage rates, means that buying a home remains extremely challenging for typical American households. Just how hard is it to afford a home at these mortgage rates? Stay tuned.

libraries and imports

labeling states with abbreviations

home value data: Zillow Home Value Index (ZHVI)

income data: median household income from FRED

combining the home value and income data

ridges and trends

”excess expensiveness”

number of states with expensiveness peaking in a given year

creating the ridgeline plot

examining the trends

conclusions

real median household income growth by state since 1984

thirties are the new twenties: US births by maternal age group

twin peaks: housing expensiveness by state, 2000-2025

diasporas and transplants: migrations between US states

syncing life stages: trends in parental age by educational attainment

getting older: children vs. total population in the US

good news about young adult parenthood rates in the US