diff --git a/openavmkit/vertical_equity_study.py b/openavmkit/vertical_equity_study.py
index b5e37286..351b4388 100644
--- a/openavmkit/vertical_equity_study.py
+++ b/openavmkit/vertical_equity_study.py
@@ -38,6 +38,26 @@ def get_vertical_equity_scores(df, sale_field: str, valuation_field: str) -> Dic
             "group_stats": None
         }
 
+    # Guard against a degenerate market proxy. When the valuation is all-zero or
+    # all-NaN (e.g. a tax-exempt model group with $0 assessed values), median_ratio
+    # collapses to 0 or NaN and market_proxy becomes non-finite; likewise a proxy
+    # with fewer distinct values than requested percentile groups makes pd.qcut
+    # raise "Bin edges must be unique". Degrade to NaN -- as we already do for
+    # fewer than 20 observations -- rather than crashing the whole modeling run.
+    market_proxy_finite = (
+        df["market_proxy"].replace([np.inf, -np.inf], np.nan).dropna()
+    )
+    if (
+        not np.isfinite(median_ratio)
+        or median_ratio == 0
+        or market_proxy_finite.nunique() < 2
+    ):
+        return {
+            "vei": np.nan,
+            "vei_significance": np.nan,
+            "group_stats": None
+        }
+
     def ci_90_lower(x):
         if len(x) <= 1: return np.nan
         # 90% CI uses the 95th percentile for a two-tailed test
@@ -47,18 +67,33 @@ def ci_90_upper(x):
         if len(x) <= 1: return np.nan
         return np.mean(x) + stats.sem(x) * stats.t.ppf(0.95, df=len(x) - 1)
 
-    #split the df into percentile_group_count groups based on the market proxy
-    df["percentile_group"] = pd.qcut(df["market_proxy"], q=percentile_group_count, labels=False)
+    # Split the df into value tiers by the market proxy. duplicates="drop"
+    # tolerates concentrated values (repeated quantile edges) instead of raising;
+    # it can yield fewer tiers than requested, so we index the actual top/bottom
+    # tiers below rather than assuming labels 0..percentile_group_count-1.
+    df["percentile_group"] = pd.qcut(
+        df["market_proxy"], q=percentile_group_count, labels=False,
+        duplicates="drop"
+    )
     grouped = df.groupby("percentile_group")
     group_stats = grouped['ratio'].agg(
         ratio='median',
         lower=ci_90_lower,
         upper=ci_90_upper
     )
-    #VEI is 100 * (median of last percentile grop - median of first percentile group)/median_ratio
-    vei_score = 100 * (group_stats[group_stats.index == (percentile_group_count - 1)].iloc[0]["ratio"] - group_stats[group_stats.index == 0].iloc[0]["ratio"]) / median_ratio
-    # VEI significance = 100 * (Lower CI for Highest PG - Upper CI for Lowest PG)/median
-    vei_significance = 100 * (group_stats[group_stats.index == (percentile_group_count -1)].iloc[0]["lower"] - group_stats[group_stats.index == 0].iloc[0]["upper"]) / median_ratio
+    if len(group_stats) < 2:
+        # need at least a bottom and top tier to measure vertical equity
+        return {
+            "vei": np.nan,
+            "vei_significance": np.nan,
+            "group_stats": None
+        }
+    bottom = group_stats.index.min()
+    top = group_stats.index.max()
+    # VEI is 100 * (top-tier median ratio - bottom-tier median ratio) / median_ratio
+    vei_score = 100 * (group_stats.loc[top, "ratio"] - group_stats.loc[bottom, "ratio"]) / median_ratio
+    # VEI significance = 100 * (Lower CI for top tier - Upper CI for bottom tier) / median_ratio
+    vei_significance = 100 * (group_stats.loc[top, "lower"] - group_stats.loc[bottom, "upper"]) / median_ratio
     return {
         "vei": vei_score,
         "vei_significance": vei_significance,
diff --git a/tests/test_vertical_equity.py b/tests/test_vertical_equity.py
new file mode 100644
index 00000000..81fbebd3
--- /dev/null
+++ b/tests/test_vertical_equity.py
@@ -0,0 +1,76 @@
+"""Tests for openavmkit.vertical_equity_study.get_vertical_equity_scores.
+
+Regression tests for a crash when the VEI "market proxy" is degenerate: with
+>= 20 observations (so the small-sample guard doesn't apply) but a market proxy
+that is all-NaN or has fewer distinct values than requested percentile groups,
+``pd.qcut`` raised ``ValueError: Bin edges must be unique``, aborting the whole
+modeling run. The function should instead degrade gracefully to a NaN result,
+exactly like it already does for < 20 observations.
+"""
+import numpy as np
+import pandas as pd
+
+from openavmkit.vertical_equity_study import get_vertical_equity_scores
+
+
+def test_all_zero_valuation_returns_nan_not_crash():
+	# Reproduces the real trigger: a tax-exempt model group whose assessed
+	# (valuation) values are all $0. Then median_ratio == 0, so
+	# market_proxy = sale*0.5 + valuation/0 = NaN for every row, and qcut on an
+	# all-NaN column raises. With >= 20 rows the small-sample guard doesn't fire.
+	rng = np.random.default_rng(0)
+	df = pd.DataFrame({
+		"sale": rng.uniform(50_000, 500_000, size=40),
+		"valuation": np.zeros(40),
+	})
+
+	result = get_vertical_equity_scores(df, "sale", "valuation")
+
+	assert np.isnan(result["vei"])
+	assert np.isnan(result["vei_significance"])
+
+
+def test_constant_market_proxy_returns_nan_not_crash():
+	# Fewer distinct market-proxy values than requested percentile groups also
+	# makes qcut raise "Bin edges must be unique". Constant sale + valuation is
+	# the extreme case (one distinct value).
+	df = pd.DataFrame({
+		"sale": np.full(30, 100_000.0),
+		"valuation": np.full(30, 90_000.0),
+	})
+
+	result = get_vertical_equity_scores(df, "sale", "valuation")
+
+	assert np.isnan(result["vei"])
+
+
+def test_concentrated_values_do_not_crash():
+	# Even with many distinct values, a heavy concentration at one value makes
+	# several quantile boundaries land on the same edge, so plain qcut(q=N) still
+	# raises "Bin edges must be unique". This is common in assessment data (capped
+	# or round-number values). Should degrade gracefully, not crash.
+	sale = np.concatenate([np.full(70, 100_000.0),
+						   np.linspace(120_000, 900_000, 30)])
+	df = pd.DataFrame({"sale": sale, "valuation": sale * 0.9})
+
+	result = get_vertical_equity_scores(df, "sale", "valuation")
+
+	# Either a finite VEI (from the tiers that could be formed) or NaN, but no raise.
+	assert np.isnan(result["vei"]) or np.isfinite(result["vei"])
+
+
+def test_healthy_data_still_returns_finite_vei():
+	# Guard must not disturb the happy path: well-spread values with a real
+	# regressive tilt should yield a finite VEI and per-group stats.
+	rng = np.random.default_rng(42)
+	true_value = rng.uniform(50_000, 500_000, size=200)
+	sale = true_value * rng.uniform(0.95, 1.05, size=200)
+	# regressive tilt: low-value assessed high, high-value assessed low
+	tilt = (true_value - true_value.min()) / (true_value.max() - true_value.min())
+	valuation = true_value * (1.15 - 0.3 * tilt)
+	df = pd.DataFrame({"sale": sale, "valuation": valuation})
+
+	result = get_vertical_equity_scores(df, "sale", "valuation")
+
+	assert np.isfinite(result["vei"])
+	assert result["group_stats"] is not None