"""A simulated residential market with known truth, for "Property Variables Most Subject to Endogeneity" (Valuation Engineer Research Labs). Every number the paper reports about bias is computed here, against a market whose true price-setting is known. Run it directly to print the results: python3 endogeneity_sim.py The market (all money in US dollars, all choices stated): L location quality the sales data do not record N(0, 1) Z size the owner started with N(0, 1) age years since construction uniform 0-80 R remodeled within the period (0/1). Owners remodel more where location is better and houses are older: P(R=1) = logistic(-1.2 + 1.0 L + 0.015 age) GLA living area, square feet: 1,700 + 300 Z + 120 L (additions are larger where location is better) C condition score, 0.00-10.00 = ten times the share of houses in the market area in worse condition, computed exactly (empirically) from a raw condition index c = -0.03 age + 1.2 R + N(0, 1); a counterfactual condition is scored against the market's actual houses. A remodel raises condition and nothing else. P price = 250,000 + 350 GLA + 90,000 L + 30,000 C - 800 age + N(0, 45,000) So the remodel has NO direct effect on price: its whole effect runs through condition (condition is the mediator). Its true effect on a given house is 30,000 x [C(R=1) - C(R=0)] for that house. Recorded data: L is never recorded. Living area from records carries an error; condition scores carry scorer noise and, in one variant, a nudge toward the known sale price (outcome-dependent measurement error). """ import numpy as np from math import erf, sqrt SEED = 20261005 N = 40_000 B0, B_GLA, B_L, B_C, B_AGE, SD_E = 250_000.0, 350.0, 90_000.0, 30_000.0, -800.0, 45_000.0 def phi(x): return 0.5 * (1.0 + np.vectorize(erf)(x / sqrt(2.0))) def market(seed=SEED, n=N): rng = np.random.default_rng(seed) L = rng.normal(0, 1, n) Z = rng.normal(0, 1, n) age = rng.uniform(0, 80, n) pr = 1.0 / (1.0 + np.exp(-(-1.2 + 1.0 * L + 0.015 * age))) R = (rng.uniform(0, 1, n) < pr).astype(float) GLA = 1_700 + 300 * Z + 120 * L noise_c = rng.normal(0, 1, n) c_raw = lambda r: -0.03 * age + 1.2 * r + noise_c # the condition score is the Journal's definition exactly: ten times the # share of the market area's houses (as they actually are) in worse # condition. A counterfactual condition is scored against the same houses. c_sorted = np.sort(c_raw(R)) C = lambda r: 10.0 * np.searchsorted(c_sorted, c_raw(r), side="left") / n C_obs_true = C(R) P = B0 + B_GLA * GLA + B_L * L + B_C * C_obs_true + B_AGE * age + rng.normal(0, SD_E, n) effect = B_C * (C(np.ones(n)) - C(np.zeros(n))) # true effect of a remodel, house by house gla_rec = GLA + rng.normal(0, 120, n) # classical error in recorded living area c_score = np.clip(C_obs_true + rng.normal(0, 0.6, n), 0, 10) # random scorer error, kept within 0-10 resid = P - (B0 + B_GLA * GLA + B_C * C_obs_true + B_AGE * age) # what the price says beyond the recorded features c_nudged = np.clip(C_obs_true + rng.normal(0, 0.6, n) + 0.8 * resid / resid.std(), 0, 10) # nudged toward the part of price the recorded features do not explain return dict(L=L, age=age, R=R, GLA=GLA, gla_rec=gla_rec, C=C_obs_true, c_score=c_score, c_nudged=c_nudged, P=P, effect=effect, C0=C(np.zeros(n)), resid=P - (B0 + B_GLA * GLA + B_L * L + B_C * C_obs_true + B_AGE * age)) def ols(y, *cols): X = np.column_stack([np.ones(len(y))] + list(cols)) b, *_ = np.linalg.lstsq(X, y, rcond=None) return b def results(m=None): m = m or market() R, P = m["R"], m["P"] out = {} out["share_remodeled"] = R.mean() out["att"] = m["effect"][R == 1].mean() # what remodels did to the houses that were remodeled out["ate"] = m["effect"].mean() # what a remodel would do to the average house out["att_unremodeled"] = m["effect"][R == 0].mean() out["raw_difference"] = P[R == 1].mean() - P[R == 0].mean() # 1. remodel flag with condition held fixed, location not recorded: the direct effect (truly 0) plus selection b = ols(P, m["GLA"], m["age"], m["C"], R); out["remodel_cond_fixed"] = b[4] # 2. remodel flag without condition (condition is the channel), location not recorded: total effect plus selection b = ols(P, m["GLA"], m["age"], R); out["remodel_total_noL"] = b[3]; out["gla_noL"] = b[1] # 3. the same with location recorded: what the sales show once the confounder is held fixed b = ols(P, m["GLA"], m["age"], R, m["L"]); out["remodel_total_withL"] = b[3]; out["gla_withL"] = b[1] # 4. living area from records (classical error), location recorded b = ols(P, m["gla_rec"], m["age"], m["C"], m["L"]); out["gla_recorded"] = b[1]; out["loc_with_gla_recorded"] = b[4] b = ols(P, m["GLA"], m["age"], m["C"], m["L"]); out["gla_true"] = b[1]; out["loc_with_gla_true"] = b[4] out["reliability_gla"] = m["GLA"].var() / m["gla_rec"].var() # 5. condition scores: classical noise vs scored knowing the price b = ols(P, m["GLA"], m["age"], m["c_score"], m["L"]); out["cond_noisy"] = b[3] b = ols(P, m["GLA"], m["age"], m["c_nudged"], m["L"]); out["cond_nudged"] = b[3] b = ols(P, m["GLA"], m["age"], m["C"], m["L"]); out["cond_true"] = b[3] b = ols(P, m["GLA"], m["age"], m["C"], R, m["L"]); out["remodel_cond_fixed_withL"] = b[4] # the remodel flag beside a condition SCORE (estimated with scorer error), location recorded: # the flag carries information about the condition the score missed b = ols(P, m["GLA"], m["age"], m["c_score"], R, m["L"]); out["remodel_beside_score_withL"] = b[4]; out["cond_beside_flag"] = b[3] b = ols(P, m["GLA"], m["age"], m["c_score"], m["L"]); out["cond_score_alone"] = b[3] # how much better is the true condition of remodeled houses than their scores say? (scores are unbiased here) out["cond_gap_remodeled"] = (m["C"][R == 1] - m["c_score"][R == 1]).mean() # age as a proxy: the age rate with condition in the model (true -800) and without it b = ols(P, m["GLA"], m["age"], m["C"], m["L"]); out["age_with_cond"] = b[2] b = ols(P, m["GLA"], m["age"], m["L"]); out["age_without_cond"] = b[2] # condition points a remodel adds, by starting condition (the rank-score mechanism) for lo, hi in ((0, 3), (3, 7), (7, 10)): k = (m["C0"] >= lo) & (m["C0"] < hi) out[f"points_from_{lo}_{hi}"] = (m["effect"][k] / B_C).mean() out[f"share_from_{lo}_{hi}"] = k.mean() out["selection_part"] = out["remodel_total_noL"] - out["att"] out["selection_part_unremodeled"] = out["remodel_total_noL"] - out["att_unremodeled"] # what a remodel does, by the condition the house starts from (the 0-10 scale's diminishing returns) for lo, hi in ((0, 3), (3, 7), (7, 10)): k = (m["C0"] >= lo) & (m["C0"] < hi) out[f"effect_from_{lo}_{hi}"] = m["effect"][k].mean() return out def simultaneity(seed=SEED, hoods=600, k=0.00003): """Neighborhood level S (e.g. a school rating) raises prices, and the neighborhood's price level raises S back (tax base, resources). P_n = 600,000 + 15,000 S_n + e_n (the effect wanted: 15,000 a point) S_n = 5 + k (P_n - 600,000) + v_n solved jointly; returns the least-squares slope of P on S against the truth.""" rng = np.random.default_rng(seed + 7) e = rng.normal(0, 60_000, hoods); v = rng.normal(0, 1.0, hoods) g = 15_000.0 # substitute: P - 600k = g (5 + k (P - 600k) + v) + e -> (P - 600k)(1 - g k) = g(5 + v) + e dP = (g * (5 + v) + e) / (1 - g * k) S = 5 + k * dP + v P = 600_000 + dP b = ols(P, S) return dict(true=g, ols=b[1], k=k, hoods=hoods) def selection(m=None, seed=SEED): """Owners sell more often when the price they can get is high for the house: the chance of a sale rises with the house's own price shock. Returns how far sold houses sit above the market's average price for their features, and the remodel/condition rates estimated on sold houses only.""" m = m or market() rng = np.random.default_rng(seed + 11) z = m["resid"] / m["resid"].std() sold = rng.uniform(0, 1, len(z)) < 1.0 / (1.0 + np.exp(-(-0.8 + 0.9 * z))) P = m["P"] b_all = ols(P, m["GLA"], m["age"], m["C"], m["L"]) b_sold = ols(P[sold], m["GLA"][sold], m["age"][sold], m["C"][sold], m["L"][sold]) level = (m["resid"][sold]).mean() return dict(share_sold=sold.mean(), level_above=level, gla_all=b_all[1], gla_sold=b_sold[1], cond_all=b_all[3], cond_sold=b_sold[3]) def theil(m=None): """The selection part of the remodel coefficient as Theil's product (entry 008, Omitted Variable Bias): gamma x delta, where gamma is what one unit of the unrecorded location quality adds to price, and delta is how much more location quality the remodeled houses have than unremodeled houses alike in their recorded features (the coefficient of the remodel flag in the auxiliary regression of location quality on the recorded features). Location quality is measured in standard deviations.""" m = m or market() R, P = m["R"], m["P"] delta = ols(m["L"], m["GLA"], m["age"], R)[3] b_short = ols(P, m["GLA"], m["age"], R)[3] b_long = ols(P, m["GLA"], m["age"], R, m["L"]) gamma_hat = b_long[4] return dict(delta=delta, gamma=B_L, gamma_hat=gamma_hat, product=gamma_hat * delta, observed_bias=b_short - b_long[3], b_short=b_short, b_long=b_long[3]) def theil_grid(t, deltas=(0.25, 0.50, 0.75, 1.00), gammas=(30_000, 60_000, 90_000, 120_000)): """The remodel gain an appraiser would conclude under assumed gamma and delta: the observed coefficient minus gamma x delta.""" return [(dl, [(g, t["b_short"] - g * dl) for g in gammas]) for dl in deltas] def spread(seeds=range(1, 11)): """Population quantities vary only slightly from one draw of 40,000 houses to another; returns the range over ten draws of the main figures.""" keys = ("att", "remodel_total_noL", "remodel_cond_fixed", "remodel_beside_score_withL", "cond_nudged", "gla_noL") vals = {k: [] for k in keys} for sd in seeds: rr = results(market(seed=sd)) for k in keys: vals[k].append(rr[k]) return {k: (min(v), max(v)) for k, v in vals.items()} def sensitivity(res, settings=(0, 20_000, 40_000, 60_000, 80_000)): """Value gained by remodeling an unremodeled house, as the appraiser's setting for the selection part varies: observed difference minus the setting.""" return [(s, res["remodel_total_noL"] - s) for s in settings] if __name__ == "__main__": r = results() for k, v in r.items(): print(f"{k:24} {v:14,.3f}" if abs(v) < 10 else f"{k:24} {v:14,.0f}") print("sensitivity:", [(s, round(v)) for s, v in sensitivity(r)]) print("simultaneity:", {k: round(v, 6) for k, v in simultaneity().items()}) print("selection:", {k: round(v, 3) for k, v in selection().items()}) # stability across seeds for sd in (1, 2, 3): rr = results(market(seed=sd)) print("seed", sd, "remodel_total_noL", round(rr["remodel_total_noL"]), "att", round(rr["att"]), "cond_nudged", round(rr["cond_nudged"]))