SPB Git
spboucher.ai

spb/wp3_uqo

Public

UQO Working Paper No. 3 — Hedonic housing price models for the US: parametric, quantile, and machine-learning approaches.

researcheconomicshousinghedonicmachine-learninguqo

R03-DATA-06/EDIT-04/ADVE-05/ADVE-12/ECON-10/THEO-08/EDIT-02/ECON-08/ADVE-17: inference v2 (fast steps) and geo-inference script

- 06_inference_v2.py: zpid/split assertions pass; v2 OLS cluster SEs: 33 of 63 slopes |t|>1.96 under state
  clustering (t(50): 33; Holm 17; BH 26), CSS effective G (state: median 14.2, min 3.7); '000' cluster = 16
  listings, nearest-ZIP3 reassignment; clustered Kennedy; lot recoding on v2 (0.0182 -> 0.0958, matches
  v1.2 Table 29); paired squared-error differences with ZIP3-cluster SEs (stored predictions)
- 09_geo_inference.py: ZIP3-cluster bootstrap of the OLS-vs-ensemble gap per design, LOSO paired tests,
  SSE decomposition, feasible level correction, 2x2 region x rounds, duplicates re-score, gap range

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Simon-Pierre Boucher committed (Nov 4, 2026) parent 52d7b12 Browse files

12 changed files +1,077 −0

12 changed files
Unified Split
added results/inference_v2/assertions.csv +4 −0
@@ -0,0 +1,4 @@
1 +,passed
2 +design_vs_sample_zpid_order,True
3 +model_data_vs_pickle_zpid_order,True
4 +ml_results_test_index_equals_seed42_split,True
added results/inference_v2/cluster_se_v2.csv +65 −0
@@ -0,0 +1,65 @@
1 +variable,coef,se_hc3,se_cluster_state,se_cluster_zip3,ratio_zip3_hc3,ratio_state_hc3,t_hc3,t_cluster_state,t_cluster_zip3,p_state_t50,p_zip3_t885,p_hc3_normal,p_holm_state,p_bh_state,p_holm_zip3,p_bh_zip3,p_holm_hc3,p_bh_hc3,G_eff_css_state,G_eff_css_zip3,n_clusters_state,n_clusters_zip3
2 +const,12.491224948743854,0.010064044077742255,0.1383233944713179,0.11084286267316142,11.013749722966997,13.744315247707963,1241.1735135748836,90.30449980269952,112.69309225237447,0.0,0.0,0.0,,,,,,,14.742894858743476,84.61237886893522,51,886
3 +ln_sqft,0.2910145212603699,0.0032281447793425463,0.0328774319743136,0.0161792092342057,5.011921812720186,10.184621267516235,90.14915412797527,8.851497935962062,17.986943431395517,8.254952277297889e-12,0.0,0.0,5.035520889151712e-10,1.7335399782325567e-10,0.0,0.0,0.0,0.0,21.080412873909594,257.89372854278923,51,886
4 +bedrooms,-0.06984165739317044,0.0022298560833404367,0.02808522300029119,0.012662345532639974,5.678548327509613,12.595083247802304,-31.32115023698934,-2.486775960171166,-5.515696694039682,0.016273719612889703,4.56126358905351e-08,0.0,0.6346750649026984,0.04100977342448205,2.189406522745685e-06,1.7959975381898197e-07,0.0,0.0,17.17117514456069,198.94941401173932,51,886
5 +bathrooms,0.2562090882976724,0.0026500747435664247,0.02347624842601837,0.015687581211757787,5.919675001560791,8.858711809171256,96.67994795983402,10.913544772925462,16.331968889228417,7.771561172376096e-15,0.0,0.0,4.89608353859694e-13,4.546363285840016e-13,0.0,0.0,0.0,0.0,14.173380814713099,142.42879365693386,51,886
6 +age,0.01221350235545477,0.013535709729939917,0.11475284332958954,0.06504842474426516,4.805689989080011,8.477785474060909,0.9023171004058599,0.10643311312448729,0.18776015566051885,0.9156646110405882,0.8511077110223486,0.3668884441604283,1.0,0.9614478415926176,1.0,0.8936630965734661,1.0,0.3852328663684497,24.319170238299606,293.67008587808834,51,886
7 +age_sq,0.044470594442945315,0.0019551074440387938,0.02093866427260368,0.008989950579174545,4.598187484061445,10.709725614541831,22.745857051763604,2.1238505887470103,4.946700657728042,0.03865239996881242,9.032448227141998e-07,0.0,1.0,0.07855165154952203,4.154926184485319e-05,3.1613568794996993e-06,0.0,0.0,25.098824545427732,397.4532214032016,51,886
8 +stories,6.9734414278304e-05,0.0006521787408528929,0.00814138935965938,0.005807423758672755,8.904650512033006,12.48337127489375,0.10692531036370206,0.00856541938945192,0.012007805384300267,0.9931999524529032,0.9904220939741246,0.9148482335305015,1.0,0.9931999524529032,1.0,0.9904220939741246,1.0,0.9296038502003483,9.636579694042908,53.506657775327284,51,886
9 +bath_per_bed,-0.027444059436691242,0.002029015687608168,0.018808759933943802,0.011785299316252374,5.808382551317309,9.26989379570339,-13.525799531418452,-1.4591105172842087,-2.3286688526310777,0.15078886599680308,0.020101007175749475,0.0,1.0,0.24358201430252804,0.5829292080967348,0.03577156571987017,0.0,0.0,10.574039482713475,84.82186338509827,51,886
10 +sqft_per_bed,-0.008940065742243996,0.0019522739180432675,0.008158415385333103,0.0063335917623245554,3.2442126608302035,4.17892966244724,-4.579309112117053,-1.0958090903677342,-1.4115317307667479,0.2784125535447177,0.15843935360323713,4.665143237492231e-06,1.0,0.3856199980829498,1.0,0.21699302776095522,6.0646862087399e-05,5.762823999255109e-06,24.554207835055568,178.23758285209254,51,886
11 +ln_lot,0.128368085196238,0.000877147580418102,0.01703129585734426,0.007551897953306306,8.609609285710635,19.41668225229112,146.34719180898833,7.537188377881588,16.998122325002683,8.696352526982309e-10,0.0,0.0,5.2178115161893857e-08,9.494743991034227e-09,0.0,0.0,0.0,0.0,17.37977584528351,282.8032746651849,51,886
12 +has_pool,0.06846104494103443,0.002607296252261303,0.05214703924529728,0.02793098670347089,10.712624880753923,20.00042733926774,26.2574860381356,1.3128462503690155,2.4510786413616747,0.19522872706864458,0.014435352347374142,0.0,1.0,0.29284309060296687,0.4763666274633467,0.0285809816869666,0.0,0.0,17.91996782502515,224.56168835932422,51,886
13 +has_spa,0.06194224195172227,0.002718639676854643,0.021506746612480983,0.01698996552718787,6.249436316196406,7.910848501028067,22.784277916294872,2.8801307360814588,3.6458132803507093,0.00583871436737593,0.00028209591748873386,0.0,0.251064717797165,0.01751614310212779,0.011848028534526822,0.0008078201273541015,0.0,0.0,12.991521495879848,175.34119498310713,51,886
14 +has_basement,-0.03657957723721726,0.002208608057829367,0.041672206518376116,0.0214857335273043,9.728178547179894,18.86808588361839,-16.562276456225504,-0.8777931454406388,-1.7025053946020294,0.3842545555416972,0.08901182530624374,0.0,1.0,0.5150646170027005,1.0,0.12744874987030355,0.0,0.0,18.222841449064315,161.46689506643904,51,886
15 +has_fireplace,0.024350781892775772,0.0018802433848401435,0.01711076631574623,0.011531290143509391,6.13287100834011,9.100293320378293,12.95086694047646,1.423126319618245,2.111713571484634,0.16090948518772197,0.03499076657170619,0.0,1.0,0.2534324391706621,0.9425370363319787,0.05801100773730237,0.0,0.0,20.521369935362646,329.65706182922094,51,886
16 +has_garage,0.11479162667680413,0.0014750242258529606,0.01604959889991434,0.010336851345120251,7.007919710025624,10.880905288612027,77.82355344735014,7.152305013517616,11.105086340531942,3.4675589066779366e-09,0.0,0.0,1.976508576806424e-07,3.120803016010143e-08,0.0,0.0,0.0,0.0,17.83111414310916,298.9306653431869,51,886
17 +on_waterfront,-0.2404703085881556,0.033099152780027404,0.4495129788578644,0.24417198192494594,7.3769858566375,13.580800144501229,-7.265149962788761,-0.5349574314831789,-0.9848398931457739,0.5950505363644241,0.32497188386153075,3.7259084706420253e-13,1.0,0.7209266113645907,1.0,0.41782099353625385,5.2162718588988355e-12,4.694644673008952e-13,5.510535873334949,62.80518907796494,51,886
18 +has_central_air,0.08348143615631584,0.0013626655549986758,0.022191812352880102,0.01394817717210821,10.23595050226522,16.28558986574051,61.26333482935728,3.761812457173262,5.985114407870542,0.0004434453476274136,3.138690640724917e-09,0.0,0.021285376686115853,0.001746066056282941,1.6007322267697077e-07,1.5210577720436135e-08,0.0,0.0,17.498804082387885,248.46294955376428,51,886
19 +has_forced_air,-0.012672739529612334,0.0015328182861712718,0.029523142761529637,0.01516491952766096,9.893488135205144,19.26069321320115,-8.26760721994435,-0.4292476458883519,-0.8356615085557911,0.6695870179105978,0.4035709046851763,2.220446049250313e-16,1.0,0.7811848542290308,1.0,0.4985287646111002,3.3306690738754696e-15,2.854859206178974e-16,28.75149520115147,345.42060323748694,51,886
20 +has_hardwood,0.12334217137791094,0.0013851055186712528,0.017653675521756323,0.00944439627051184,6.818539196617998,12.74536508864082,89.04893505603418,6.986770048305494,13.059825937526698,6.293867205187098e-09,0.0,0.0,3.524565634904775e-07,4.95642042408484e-08,0.0,0.0,0.0,0.0,20.821706230738403,335.4267238273962,51,886
21 +parking_spaces,0.021179080148683513,0.0007340355349200892,0.005519964544873688,0.004314311335556711,5.877523812285721,7.520023598686703,28.852935779177468,3.8368145259831823,4.909029159331776,0.00035078728722570673,1.0896243947122741e-06,0.0,0.01718857707405963,0.0014733066063479683,4.9033097762052336e-05,3.612965098256488e-06,0.0,0.0,21.104708218715135,208.01235284090936,51,886
22 +luxury_score,0.004193343993685188,0.0019022725698461488,0.008848805270257743,0.008549381768231569,4.4942990314595255,4.6517020801984295,2.2043865112476158,0.47388815389346306,0.49048505580451784,0.637642516816556,0.6239122843888503,0.027497173481998427,1.0,0.7579524256498684,1.0,0.7278976651203254,0.19326922360586174,0.030391612795893003,18.02107906095303,268.4680953559843,51,886
23 +walk_score,0.03662780227257727,0.0011750988502931105,0.01918261650586889,0.012397583396033071,10.550247234894904,16.324257743154188,31.169975413933066,1.9094268115807376,2.9544308033690894,0.06195337445771276,0.003215453008393432,0.0,1.0,0.11479595855399718,0.12218721431895041,0.007791289981876392,0.0,0.0,18.86209219667965,316.4946705014394,51,886
24 +bike_score,0.07865281426594363,0.0009971544190365488,0.027131915800962563,0.012152014141512105,12.186692361302864,27.209342187118256,78.87726591227266,2.898903816558109,6.472409705092448,0.00554820169142034,1.5958456778264463e-10,0.0,0.24412087442249497,0.01747683532797407,8.457982092480165e-09,9.139843427551465e-10,0.0,0.0,20.289265360271763,276.623578644146,51,886
25 +transit_score,0.05398709618488032,0.0008091091712144606,0.012421369463718265,0.01809145703790151,22.359723114677475,15.351907882930035,66.72411845715023,4.346307896449897,2.984120962273941,6.786283914395064e-05,0.002922111645333736,0.0,0.0035967304746293838,0.0003886689878244446,0.1139623541680157,0.007363721346241014,0.0,0.0,14.229812245284245,83.27509061842643,51,886
26 +avg_school_rating,-0.047799518130700575,0.0007052258546274312,0.023743534725432027,0.013389621998098181,18.986289158629727,33.66798674443901,-67.77902117039218,-2.0131593161443586,-3.569893021437787,0.04949729596693664,0.00037634780399264933,0.0,1.0,0.09449483775506086,0.015430259963698623,0.0010308657239798656,0.0,0.0,15.266131686280634,110.9291394233956,51,886
27 +school_count,-0.000755804398075912,0.0005117978849218512,0.0034809827641059743,0.002510806710253446,4.905855972102801,6.801479386022671,-1.4767634262328364,-0.2171238553288345,-0.3010205425170381,0.8289953823790519,0.7634695367912179,0.13973906359753863,1.0,0.9162580542084259,1.0,0.8292858761697712,0.6986953179876931,0.14921289841771074,8.190317769272221,44.51533969050612,51,886
28 +nearest_school_distance,-0.014084634096873461,0.0007335863904315645,0.012938684630471637,0.005430923137351372,7.403249580675007,17.637574523240385,-19.19969383372496,-1.0885676944086744,-2.5934143681772763,0.28156380812405857,0.009659761825693725,0.0,1.0,0.3856199980829498,0.33809166389928036,0.020984999828231196,0.0,0.0,13.91941393201027,101.25645034902385,51,886
29 +property_tax_rate,-0.018914530233367303,0.000777317956299977,0.026918130276302472,0.014635447354455232,18.828134916784588,34.62949756677745,-24.333067414781222,-0.7026687975434467,-1.2923780035742782,0.48552092517498013,0.1965636036937275,0.0,1.0,0.6372462142921614,1.0,0.26347887303627304,0.0,0.0,14.80823639322118,189.78537010701427,51,886
30 +is_condo,-0.09331344936292597,0.003737929726529345,0.054374020080158014,0.040180853579487565,10.749494110151543,14.546560277537402,-24.963938915343707,-1.7161403410188096,-2.3223361638729036,0.09232696683214003,0.020440894697068668,0.0,1.0,0.16618854029785207,0.5829292080967348,0.03577156571987017,0.0,0.0,7.808686326022657,106.35370988349715,51,886
31 +has_hoa,-0.06580598259352977,0.0023744807864962175,0.02872201278520274,0.02246088032821604,9.459280721895965,12.096123476149465,-27.71384083955172,-2.291134088883641,-2.9298042477374446,0.026205928787518484,0.003478894600409621,0.0,0.9434134363506654,0.05896333977191659,0.12871910021515598,0.008117420734289116,0.0,0.0,13.431641567799643,131.20947361770786,51,886
32 +ln_hoa,0.057869838061030084,0.0013934965179893742,0.0193044246353199,0.012635782200760275,9.067681216019103,13.853227751995792,41.52851285522291,2.997749953922473,4.579838204044712,0.004228020187029946,5.319860956198852e-06,0.0,0.19026090841634757,0.014019224830678243,0.00023407388207274948,1.6757562012026384e-05,0.0,0.0,11.087480568317162,112.17215606286047,51,886
33 +tag_new_construction,0.07107576583056167,0.0029117712367649345,0.030205497245712504,0.02278278972506957,7.82437488131174,10.37358184747773,24.40980422264525,2.353073854483573,3.1197130240969484,0.022593443011505343,0.0018689191259826377,0.0,0.8359573914256977,0.05271803369351247,0.0747567650393055,0.004905912705704424,0.0,0.0,13.11778060251481,204.75897125079314,51,886
34 +tag_foreclosure,-0.3895224827615263,0.009701797276632012,0.036313094915301145,0.022168248586423734,2.2849630799665053,3.742924520054213,-40.14951783209692,-10.726777314631871,-17.571188866948987,1.4432899320127035e-14,0.0,0.0,8.948397578478762e-13,4.546363285840016e-13,0.0,0.0,0.0,0.0,14.907021133619281,163.98534085042598,51,886
35 +sqft_x_age,-0.1583055650153391,0.013557565877695544,0.1079391947593836,0.06464022434229162,4.767834058518983,7.961546765335035,-11.676547725707765,-1.466617991437044,-2.449025612551376,0.14874174429546105,0.014517324031475098,0.0,1.0,0.24358201430252804,0.4763666274633467,0.0285809816869666,0.0,0.0,25.06952124574614,303.7508126748135,51,886
36 +pool_x_south,0.004673129995258252,0.0011656355783042935,0.030280776327294562,0.015360411081263122,13.177712972358528,25.977910155544045,4.009083183662325,0.154326624415041,0.30423209187146116,0.8779734599515412,0.7610225961099109,6.095494332503115e-05,1.0,0.9374970843550355,1.0,0.8292858761697712,0.0007314593199003738,7.384925825917235e-05,17.942790592842513,251.87353761845287,51,886
37 +waterfront_x_sqft,0.17142561756808375,0.009319182779343314,0.12436199177497348,0.06815329712442604,7.313226785882241,13.344731476952186,18.39491955754552,1.3784405920280647,2.5152945609530177,0.17420414973441467,0.012069431445146828,0.0,1.0,0.26767954715288106,0.41036066913499214,0.02534580603480834,0.0,0.0,5.819036279992954,68.04236420332771,51,886
38 +basement_x_north,-0.02291213287340406,0.0011737579744655828,0.018914924249803747,0.012133956179348968,10.337698608500284,16.114841953185277,-19.520321371052713,-1.2113256479809473,-1.8882656682408903,0.23146471739851937,0.05931715290533002,0.0,1.0,0.3391227254908539,1.0,0.09114586909843393,0.0,0.0,28.417657810141886,336.3739303819459,51,886
39 +condo_x_walkscore,0.037280704479362194,0.0013630021654134881,0.017119127845219746,0.01558509183297419,11.43438523316374,12.559868413728008,27.351904072766096,2.1777221840055514,2.3920747390455195,0.03417141108346322,0.016961210121812842,0.0,1.0,0.07175996327527277,0.5257975137761981,0.031993839941966785,0.0,0.0,9.45155781981122,93.09507504199935,51,886
40 +age_x_luxury,0.036682706662375185,0.0014405003826735725,0.01521064902082495,0.00827906339835743,5.747352446371077,10.559281485641778,25.465252979865223,2.4116463809106876,4.430779775119617,0.0195943833197747,1.0565875494750188e-05,0.0,0.7445865661514386,0.04747869804406946,0.00045433264627425807,3.169762648425056e-05,0.0,0.0,23.51502079033498,380.34247949816586,51,886
41 +roof_cat_Metal,0.018216080283262076,0.008079529277938074,0.032030972793119325,0.04349388062002295,5.383219631221232,3.964460266340386,2.254596729168718,0.5687020622481728,0.4188193838669819,0.5721040357970757,0.6754497603961143,0.024158652950732717,1.0,0.70671675010227,1.0,0.7698064415575635,0.19326922360586174,0.027178484569574307,11.597920195310369,49.49503921869862,51,886
42 +roof_cat_Other,0.025633184210792352,0.008020183898511508,0.04051829676988926,0.03934087130984624,4.905233072915901,5.052040861233765,3.1960843460895823,0.6326323230309469,0.6515662555846056,0.5298564588285082,0.5148501808159256,0.001393063818063256,1.0,0.6812440184937962,1.0,0.6237607959885252,0.01393063818063256,0.0016252411210737987,11.2209981565301,49.05588359006198,51,886
43 +roof_cat_Shingle_Comp,-0.08722454365265524,0.007672273400325394,0.05841358500107941,0.04194027562623559,5.466473030595708,7.613595339107961,-11.368800237092163,-1.4932235994596708,-2.079732246635314,0.14166259366329026,0.037837438945826474,0.0,1.0,0.2412092811023591,0.9425370363319787,0.059593966339676695,0.0,0.0,11.148440948664517,43.45468319600791,51,886
44 +roof_cat_Slate,-0.001479329275050758,0.017354319316023003,0.07048985061618178,0.05103160297577027,2.940570704415557,4.061804403420161,-0.08524271382311802,-0.020986415237361285,-0.028988493184373243,0.983340019856497,0.9768803034149873,0.9320684334169165,1.0,0.9931999524529032,1.0,0.9904220939741246,1.0,0.9320684334169165,17.627100473408593,103.77458879111367,51,886
45 +roof_cat_Tile,-0.10655882776595488,0.007946053630680953,0.06567413623045368,0.04466982958354879,5.621637061581311,8.265000374132335,-13.410282980536982,-1.622538702176986,-2.3854764784058826,0.11097738655256939,0.017266516794077313,0.0,1.0,0.19421042646699643,0.5257975137761981,0.031993839941966785,0.0,0.0,10.457749432549283,45.06254488521472,51,886
46 +roof_cat_Unknown,-0.02498087593032834,0.007709418740984999,0.07128071539237024,0.04770444842551859,6.187813897293066,9.245926027266618,-3.240306016525526,-0.350457705044333,-0.5236592551600553,0.7274674120507936,0.6006466607249532,0.0011940148283351437,1.0,0.8332808538036363,1.0,0.7139762193523028,0.013134163111686581,0.001419300645002152,11.276093233914358,44.22757492689323,51,886
47 +construction_cat_Concrete_Block,0.2390621238483508,0.0027062872326422105,0.031745828637042166,0.03682500870417783,13.607206308335842,11.730398848332142,88.33582812824673,7.530505080891308,6.491841611465226,8.907317106121582e-10,1.4111689594642485e-10,0.0,5.255317092611733e-08,9.494743991034227e-09,7.620312381106942e-09,8.890364444624765e-10,0.0,0.0,4.256006089072552,61.92572800728507,51,886
48 +construction_cat_Other,0.06554233625298557,0.002245038521643814,0.02125220544078923,0.01291572716024121,5.7530091513905655,9.466298789932743,29.194303626022222,3.084025158499106,5.074614494393015,0.003321785503759944,4.732546519115033e-07,0.0,0.15280213317295743,0.011626249263159805,2.2242968639840655e-05,1.7538260629661594e-06,0.0,0.0,19.647810299396067,257.6525888779379,51,886
49 +construction_cat_Steel,0.005430789119110303,0.01301629395330231,0.06852248782141047,0.038287788169592996,2.9415276196861826,5.264362349778211,0.41723006092164044,0.0792555742176863,0.1418412861838614,0.9371456981913953,0.8872376539883065,0.6765101393150521,1.0,0.9678717866566869,1.0,0.9163274131354641,1.0,0.6986907996204637,7.703413398151516,41.958414112151495,51,886
50 +construction_cat_Stone,0.13309343926465478,0.004418217902119623,0.017683776131302405,0.014723890394726274,3.3325405674678343,4.002468081716541,30.123783437843503,7.52630197738499,9.0392848422945,9.042613324794502e-10,0.0,0.0,5.255317092611733e-08,9.494743991034227e-09,0.0,0.0,0.0,0.0,11.226577002888835,120.88907199298262,51,886
51 +construction_cat_Stucco,0.15173158300246475,0.0020508506260864742,0.032983581236229124,0.02606935999504513,12.711486474659473,16.08287840015432,73.98470716124548,4.600215541052376,5.820303338144995,2.9089775380075622e-05,8.20905499132607e-09,0.0,0.0015708478705240836,0.0001832655848944764,4.0224369457497744e-07,3.4478030963569496e-08,0.0,0.0,8.369079399059057,130.5060305566327,51,886
52 +construction_cat_Unknown,0.17079441684902644,0.001999999449781519,0.044336863753570314,0.026759383209593708,13.379695285674662,22.16843797552729,85.39723191808083,3.852198878980764,6.3825991619939,0.00033423580879232695,2.805715659803809e-10,0.0,0.017046026248408674,0.0014733066063479683,1.4589721430979807e-08,1.4730007213969998e-09,0.0,0.0,14.850779836585065,200.56217009022106,51,886
53 +construction_cat_Vinyl,0.004651255255503686,0.0017632186688960332,0.020802870961154317,0.011434497065796292,6.485013610339966,11.798236559155297,2.637934442025774,0.2235871800670726,0.4067739253190997,0.8239895719534498,0.684272392495612,0.008341270561836645,1.0,0.9162580542084259,1.0,0.7698064415575635,0.0750714350565298,0.009554546279921975,25.03031930949464,335.6784841298451,51,886
54 +construction_cat_Wood_Frame,0.09008496696394888,0.002100159011884044,0.032346199632941275,0.015193432991140286,7.234420301113447,15.401785983778263,42.894355357946935,2.7850247629154743,5.92920421714302,0.007539783862656746,4.360506622447247e-09,0.0,0.31667092223158333,0.020806762911823193,2.1802533112236233e-07,1.962227980101261e-08,0.0,0.0,24.197366634706995,338.0855527564923,51,886
55 +foundation_cat_Concrete,0.20961928119498374,0.005002852887764971,0.04518270560418669,0.07990891509398913,15.972669372192655,9.031388013564417,41.89994906858661,4.639369829494234,2.6232277205669585,2.5490642459269353e-05,0.00885996297263203,0.0,0.0014019853352598144,0.00017843449721488547,0.31895866701475306,0.019934916688422066,0.0,0.0,4.222787367216655,25.57745619019931,51,886
56 +foundation_cat_Crawlspace,0.1730141087771627,0.005596998631607597,0.062185442734761204,0.08189368995768126,14.63171520093071,11.110498112968093,30.91194409091159,2.782228463261372,2.1126671525799865,0.0075961197932052915,0.03490877912340662,0.0,0.31667092223158333,0.020806762911823193,0.9425370363319787,0.05801100773730237,0.0,0.0,5.942448041597181,37.15132914004164,51,886
57 +foundation_cat_Other,0.1420168557497226,0.005679964997896968,0.054319395695307596,0.08176791319712322,14.395848077831138,9.563332822547249,25.003121639359563,2.6144778293619857,1.7368286678339624,0.011778963940053844,0.08276536290315772,0.0,0.47115855760215375,0.03091978034264134,1.0,0.12126088053253341,0.0,0.0,6.315538118068315,38.07586039607309,51,886
58 +foundation_cat_Pier_Raised,0.15687380700563638,0.006255611885023543,0.07488609703107221,0.08650916880832156,13.829049883262696,11.97102672088014,25.07729218003524,2.0948321948271027,1.8133778091571058,0.04127309683048397,0.0701121916115568,0.0,1.0,0.08125640938501531,1.0,0.1051682874173352,0.0,0.0,6.8404516250198135,45.51784567567278,51,886
59 +foundation_cat_Slab,0.12780073865115896,0.005039709078499911,0.05713150096585598,0.07981374203303475,15.836974077240113,11.336269629051959,25.35875318604687,2.2369574838850754,1.6012372731285107,0.02977814528251921,0.10968124089834719,0.0,1.0,0.06469045354478312,1.0,0.15355373725768606,0.0,0.0,4.302769522447288,25.935817352264262,51,886
60 +foundation_cat_Stone,0.02066416003543145,0.00959212726767504,0.12966557486527283,0.08938436025355444,9.318512751053145,13.5179164378103,2.1542833470390423,0.1593650439363127,0.23118317317273324,0.8740232957518832,0.817225938034436,0.031217958024114756,1.0,0.9374970843550355,1.0,0.8726310863757537,0.19326922360586174,0.033909161302055686,11.743390086617975,110.53874045204155,51,886
61 +foundation_cat_Unknown,0.1645492571308587,0.004903396729553319,0.04725545744805438,0.07895862665340403,16.102842786004157,9.637290240710174,33.55821815092015,3.482121769993226,2.083993403952711,0.0010425908170186293,0.03744710804794349,0.0,0.04900176839987558,0.003863718910127862,0.9425370363319787,0.059593966339676695,0.0,0.0,3.663564544662053,22.031024576681364,51,886
62 +region_Northeast,0.4765442256871757,0.003007095699199933,0.12383183900188474,0.05357587517963916,17.816484920614112,41.17987965425623,158.4732490602028,3.8483174402337883,8.89475391842559,0.00033833928627258203,0.0,0.0,0.017046026248408674,0.0014733066063479683,0.0,0.0,0.0,0.0,12.508353076597146,231.9655557431667,51,886
63 +region_South,0.03509993423093495,0.002672271598294191,0.061614473439691904,0.03385698865833608,12.669740860153674,23.056965272176182,13.134867823068781,0.569670278287628,1.0367116398077194,0.5714520837525607,0.3001534413119713,0.0,1.0,0.70671675010227,1.0,0.39395139172196236,0.0,0.0,19.370165805003577,270.1900666119216,51,886
64 +region_West,0.41948338327278023,0.0030348961805360514,0.10424876957313231,0.047435232182475966,15.629935708079977,34.35002826182968,138.2200109391179,4.02386891462067,8.843287235510788,0.00019380770165189354,0.0,0.0,0.010078000485898464,0.001017490433672441,0.0,0.0,0.0,0.0,21.09398621373925,304.641878093168,51,886
65 +lot_missing,-0.023018541083695682,0.0023514230399513717,0.020192316483941774,0.026876878999291056,11.430048333560123,8.587275084435406,-9.789196028364048,-1.1399653478094751,-0.8564439749236826,0.25973226425867546,0.39198408126009454,0.0,1.0,0.3718893783703762,1.0,0.4938999423877191,0.0,0.0,12.280714701251377,171.76161296859078,51,886
added results/inference_v2/kennedy_clustered_v2.csv +40 −0
@@ -0,0 +1,40 @@
1 +variable,se,coef,se_value,naive_pct,kennedy_pct,t
2 +on_waterfront,hc3,-0.2404703085880398,0.03309915278002476,-21.374200978856006,-21.417258585756105,-7.265149962785842
3 +on_waterfront,cluster_zip3,-0.2404703085880398,0.2441719819246039,-21.374200978856006,-23.683444100288476,-0.9848398931466793
4 +on_waterfront,cluster_state,-0.2404703085880398,0.44951297885765096,-21.374200978856006,-28.9297433316168,-0.5349574314831752
5 +has_pool,hc3,0.06846104494106407,0.002607296252260914,7.085890876655343,7.085526892665639,26.257486038150883
6 +has_pool,cluster_zip3,0.06846104494106407,0.02793098670346923,7.085890876655343,7.044128027971941,2.451078641362882
7 +has_pool,cluster_state,0.06846104494106407,0.05214703924533937,7.085890876655343,6.9403897496938205,1.3128462503685243
8 +has_garage,hc3,0.11479162667678366,0.0014750242258530118,12.163969342939108,12.163847325629561,77.82355344733355
9 +has_garage,cluster_zip3,0.11479162667678366,0.010336851345121715,12.163969342939108,12.157977115144924,11.10508634052839
10 +has_garage,cluster_state,0.11479162667678366,0.016049598899915894,12.163969342939108,12.149524135800682,7.152305013515647
11 +has_central_air,hc3,0.08348143615631719,0.0013626655549987552,8.706503500018314,8.706402573826644,61.2633348293547
12 +has_central_air,cluster_zip3,0.08348143615631719,0.013948177172114098,8.706503500018314,8.695929499707455,5.985114407868112
13 +has_central_air,cluster_state,0.08348143615631719,0.022191812352839208,8.706503500018314,8.679739094248816,3.7618124571802545
14 +has_hardwood,hc3,0.12334217137788535,0.0013851055186713268,13.127144343826735,13.12703582568715,89.04893505601096
15 +has_hardwood,cluster_zip3,0.12334217137788535,0.009444396270530644,13.127144343826735,13.122099176825941,13.059825937497987
16 +has_hardwood,cluster_state,0.12334217137788535,0.017653675521770232,13.127144343826735,13.109517552151084,6.986770048298539
17 +is_condo,hc3,-0.09331344936292371,0.003737929726530545,-8.90920679680941,-8.909843160271425,-24.96393891533509
18 +is_condo,cluster_zip3,-0.09331344936292371,0.04018085357954053,-8.90920679680941,-8.982710213085909,-2.322336163869786
19 +is_condo,cluster_state,-0.09331344936292371,0.05437402008008411,-8.90920679680941,-9.043763833004403,-1.7161403410211005
20 +tag_foreclosure,hc3,-0.38952248276153817,0.009701797276631232,-32.261974196838224,-32.26516203827369,-40.149517832101374
21 +tag_foreclosure,cluster_zip3,-0.38952248276153817,0.022168248586454463,-32.261974196838224,-32.2786164433149,-17.571188866925166
22 +tag_foreclosure,cluster_state,-0.38952248276153817,0.03631309491569104,-32.261974196838224,-32.306620541476995,-10.726777314517024
23 +tag_new_construction,hc3,0.07107576583051099,0.002911771236765836,7.366256972305973,7.365801825604157,24.409804222620284
24 +tag_new_construction,cluster_zip3,0.07107576583051099,0.022782789725085007,7.366256972305973,7.338396064292785,3.11971302409261
25 +tag_new_construction,cluster_state,0.07107576583051099,0.03020549724569273,7.366256972305973,7.317289155640827,2.3530738544834358
26 +has_hoa,hc3,-0.06580598259351089,0.0023744807864974193,-6.368749240554539,-6.369013194122319,-27.713840839529748
27 +has_hoa,cluster_zip3,-0.06580598259351089,0.022460880328242394,-6.368749240554539,-6.3923643304857,-2.9298042477331667
28 +has_hoa,cluster_state,-0.06580598259351089,0.02872201278522879,-6.368749240554539,-6.407362014848806,-2.2911340888809058
29 +has_basement,hc3,-0.03657957723720858,0.002208608057829552,-3.591862808831825,-3.5920979455599467,-16.56227645622019
30 +has_basement,cluster_zip3,-0.03657957723720858,0.021485733527307798,-3.591862808831825,-3.6141130101875807,-1.7025053946013482
31 +has_basement,cluster_state,-0.03657957723720858,0.04167220651835502,-3.591862808831825,-3.6755363514602513,-0.8777931454408748
32 +has_spa,hc3,0.06194224195172405,0.0027186396768545877,6.390089411421385,6.3896962474824015,22.78427791629599
33 +has_spa,cluster_zip3,0.06194224195172405,0.016989965527204557,6.390089411421385,6.3747352948613445,3.645813280347233
34 +has_spa,cluster_state,0.06194224195172405,0.021506746612520354,6.390089411421385,6.36548741243439,2.8801307360762687
35 +has_fireplace,hc3,0.02435078189277515,0.0018802433848400673,2.464968341285667,2.4647872184664132,12.950866940476654
36 +has_fireplace,cluster_zip3,0.02435078189277515,0.01153129014352735,2.464968341285667,2.4581561509005923,2.111713571481291
37 +has_fireplace,cluster_state,0.02435078189277515,0.017110766315633092,2.464968341285667,2.4499696782881664,1.4231263196276187
38 +lot_missing,hc3,-0.023018541083680823,0.002351423039951501,-2.27556355642361,-2.2758337245538107,-9.789196028357189
39 +lot_missing,cluster_zip3,-0.023018541083680823,0.026876878999282296,-2.27556355642361,-2.31085361860347,-0.856443974923409
40 +lot_missing,cluster_state,-0.023018541083680823,0.020192316483956356,-2.27556355642361,-2.2954841006944116,-1.139965347807916
added results/inference_v2/lot_recoding_v2.csv +11 −0
@@ -0,0 +1,11 @@
1 +sample,specification,N,R2,beta_ln_lot,se_ln_lot,beta_ln_sqft,beta_lot_missing
2 +"all (788,842)","v1.1 coding: ln(lot) = 0 when missing, zero, 1 sq ft or negative",788842,0.6342771060419653,0.018228359685209458,0.00025052415374509794,0.6280515162728741,
3 +"all (788,842)",v2: ln(state-median lot) when missing/zero/<= 1,788842,0.6482128957511863,0.09580171656077102,0.0006800775296213197,0.5849975385535291,
4 +"all (788,842)",v2 + lot-missing indicator (the design used throughout),788842,0.6482563185065394,0.09590009318528735,0.0006805181853858852,0.5865087924645662,-0.021332999595101404
5 +"all (788,842)",Complete cases (raw lot > 1 sq ft),658571,0.660083194475785,0.09740869847290046,0.0007011909144493823,0.545118686867461,
6 +"all (788,842)","v2 + indicator, lot winsorised at p99.5",788842,0.6478467301733439,0.09705411277895361,0.0006600607717048909,0.5854917282575403,-0.02209064164386547
7 +"trimmed (783,872)","v1.1 coding: ln(lot) = 0 when missing, zero, 1 sq ft or negative",783872,0.6355094592975552,0.017589647583908158,0.00024992370757992796,0.6192778822625188,
8 +"trimmed (783,872)",v2: ln(state-median lot) when missing/zero/<= 1,783872,0.6486347206805829,0.09668354003930296,0.000661218056540847,0.5851520876966476,
9 +"trimmed (783,872)",v2 + lot-missing indicator (the design used throughout),783872,0.6486853164598687,0.09681385576466428,0.0006615354525631342,0.5867712171353324,-0.023018541083680823
10 +"trimmed (783,872)",Complete cases (raw lot > 1 sq ft),655220,0.6612844635289765,0.09818146519520987,0.0006836538611093285,0.5363169456576322,
11 +"trimmed (783,872)","v2 + indicator, lot winsorised at p99.5",783872,0.6480153086487406,0.09609756392239553,0.0006595274427533722,0.5864874716437819,-0.023284590576429752
added results/inference_v2/ols_v2_unstandardized_cluster_se.csv +65 −0
@@ -0,0 +1,65 @@
1 +variable,coef,se_hc3,se_cluster_zip3,se_cluster_state
2 +const,7.097385173820015,0.04713869307832596,0.25978597052235713,0.40909620176511574
3 +ln_sqft,0.5867712171353324,0.006508893209382003,0.03262206385566403,0.06629061220461152
4 +bedrooms,-0.06429849088151064,0.002052877700688479,0.011657365233845963,0.02585616553730326
5 +bathrooms,0.22127205088268748,0.002288706764451514,0.013548400219604724,0.020274993642016063
6 +age,0.0003993303401917053,0.0004425609799430484,0.002126810870966191,0.003751937047033228
7 +age_sq,1.1798857925959607e-05,5.187255815152157e-07,2.38519747658669e-06,5.555408647269214e-06
8 +stories,0.00011189635296935936,0.0010464907942168474,0.009318634786547892,0.013063733119953604
9 +bath_per_bed,-0.09629919767750475,0.0071196676731510025,0.04135375348428236,0.06599856319056782
10 +sqft_per_bed,-3.5798930813020545e-05,7.817539706607854e-06,2.536176129495805e-05,3.2668948563922504e-05
11 +ln_lot,0.09681385576466428,0.0006615354525631343,0.0056955617752133285,0.012844823681033703
12 +has_pool,0.06846104494106407,0.002607296252260914,0.02793098670346923,0.05214703924533937
13 +has_spa,0.06194224195172405,0.0027186396768545877,0.016989965527204557,0.021506746612520354
14 +has_basement,-0.03657957723720858,0.002208608057829552,0.021485733527307798,0.04167220651835502
15 +has_fireplace,0.02435078189277515,0.0018802433848400673,0.01153129014352735,0.017110766315633092
16 +has_garage,0.11479162667678366,0.0014750242258530118,0.010336851345121715,0.016049598899915894
17 +on_waterfront,-0.2404703085880398,0.03309915278002476,0.2441719819246039,0.44951297885765096
18 +has_central_air,0.08348143615631719,0.0013626655549987552,0.013948177172114098,0.022191812352839208
19 +has_forced_air,-0.012672739529595476,0.0015328182861712425,0.015164919527650932,0.029523142761524648
20 +has_hardwood,0.12334217137788535,0.0013851055186713268,0.009444396270530644,0.017653675521770232
21 +parking_spaces,0.015274264289014274,0.0005293833669444296,0.003111463345039875,0.0039809754121526465
22 +luxury_score,0.0032520835165942893,0.0014752782690289515,0.006630341695593732,0.006862554992930626
23 +walk_score,0.0013836486341733901,4.439043072039199e-05,0.00046833001896442675,0.0007246408324096627
24 +bike_score,0.004359040419130759,5.526358411034659e-05,0.0006734802983364233,0.001503685770547418
25 +transit_score,0.0029608941661074745,4.4375170996232346e-05,0.0009922165366431797,0.0006812435374233082
26 +avg_school_rating,-0.05289868104100902,0.0007804580256189923,0.014818001750603289,0.026276450461196633
27 +school_count,-0.008802206509876473,0.005960471632475837,0.029241215356429567,0.040540024939686434
28 +nearest_school_distance,-0.023067406252163858,0.001201446567428564,0.008894608796531407,0.021190603368743635
29 +property_tax_rate,-0.03588590697177712,0.0014747794168354436,0.02776734583269444,0.05107087022681702
30 +is_condo,-0.09331344936292371,0.003737929726530545,0.04018085357954053,0.05437402008008411
31 +has_hoa,-0.06580598259351089,0.0023744807864974193,0.022460880328242394,0.02872201278522879
32 +ln_hoa,0.015470542525758964,0.00037252820922590186,0.0033779670452368463,0.0051607181264401
33 +tag_new_construction,0.07107576583051099,0.002911771236765836,0.022782789725085007,0.03020549724569273
34 +tag_foreclosure,-0.38952248276153817,0.009701797276631232,0.022168248586454463,0.03631309491569104
35 +sqft_x_age,-0.0007009127976678083,6.0027399718939595e-05,0.0002862006808364781,0.00047791095005187057
36 +pool_x_south,0.01041965273957369,0.0025990113605197343,0.03424902572083211,0.06751688361683991
37 +waterfront_x_sqft,0.0802433137978051,0.004362254129287554,0.03190215374509432,0.05821310998955928
38 +basement_x_north,-0.06297566715394295,0.003226159342200936,0.033351062942684986,0.051989047915319046
39 +condo_x_walkscore,0.0017934635921976702,6.556997229252562e-05,0.00074975232292117,0.0008235502238869865
40 +age_x_luxury,0.0005962715057854183,2.3415102385084608e-05,0.0001345748459758663,0.0002472466570972705
41 +roof_cat_Metal,0.01821608028351448,0.00807952927793866,0.0434938806200699,0.03203097279258099
42 +roof_cat_Other,0.0256331842110276,0.008020183898512152,0.03934087130992595,0.04051829676947301
43 +roof_cat_Shingle_Comp,-0.08722454365244332,0.007672273400325875,0.04194027562626844,0.05841358500067677
44 +roof_cat_Slate,-0.0014793292747906361,0.01735431931602481,0.05103160297588372,0.07048985061580963
45 +roof_cat_Tile,-0.10655882776575358,0.007946053630681673,0.0446698295835965,0.06567413623023849
46 +roof_cat_Unknown,-0.02498087593011189,0.007709418740985464,0.04770444842552367,0.07128071539186083
47 +construction_cat_Concrete_Block,0.23906212384833983,0.0027062872326424656,0.03682500870428442,0.03174582863705286
48 +construction_cat_Other,0.06554233625299936,0.0022450385216437408,0.012915727160244209,0.021252205440810714
49 +construction_cat_Steel,0.00543078911915027,0.013016293953302638,0.0382877881696356,0.06852248782173655
50 +construction_cat_Stone,0.13309343926463582,0.004418217902119798,0.014723890394724612,0.01768377613128725
51 +construction_cat_Stucco,0.1517315830024574,0.0020508506260866633,0.02606935999505983,0.03298358123639612
52 +construction_cat_Unknown,0.1707944168490404,0.0019999994497814505,0.026759383209622473,0.04433686375365333
53 +construction_cat_Vinyl,0.004651255255493689,0.0017632186688960516,0.011434497065813629,0.020802870961178527
54 +construction_cat_Wood_Frame,0.09008496696394619,0.0021001590118839312,0.015193432991165264,0.03234619963305997
55 +foundation_cat_Concrete,0.20961928119512518,0.005002852887764494,0.0799089150938967,0.04518270560414532
56 +foundation_cat_Crawlspace,0.17301410877727674,0.005596998631607147,0.08189368995758226,0.06218544273471564
57 +foundation_cat_Other,0.1420168557498723,0.005679964997896459,0.08176791319701499,0.05431939569520319
58 +foundation_cat_Pier_Raised,0.15687380700578496,0.006255611885023287,0.0865091688082324,0.07488609703095109
59 +foundation_cat_Slab,0.127800738651305,0.005039709078499568,0.07981374203293969,0.057131500965805576
60 +foundation_cat_Stone,0.02066416003562676,0.009592127267675719,0.08938436025354381,0.1296655748654012
61 +foundation_cat_Unknown,0.16454925713098043,0.0049033967295528304,0.07895862665330432,0.04725545744780093
62 +region_Northeast,0.47654422568710914,0.003007095699199743,0.05357587517968932,0.12383183900187825
63 +region_South,0.03509993423092984,0.002672271598293785,0.03385698865837279,0.06161447343962
64 +region_West,0.4194833832728059,0.003034896180535963,0.0474352321825213,0.1042487695733401
65 +lot_missing,-0.023018541083680823,0.002351423039951501,0.026876878999282296,0.020192316483956356
added results/inference_v2/paired_loss_cluster_se.csv +4 −0
@@ -0,0 +1,4 @@
1 +predictions,quantity,mean,se_iid,se_hc0,se_cluster_zip3,t_cluster_zip3,n_clusters
2 +v1.2 stored,OLS minus XGBoost squared-error difference,0.13230822436471384,0.0009393288059222971,0.0009393258290057469,0.006103686800595308,21.676771545979307,885
3 +v1.2 stored,LightGBM minus XGBoost squared-error difference,0.015807840336707984,0.0002180070873428965,0.00021800639643583365,0.0005837125659902828,27.08154879258012,885
4 +v1.2 stored,OLS minus LightGBM squared-error difference,0.11650038402800585,0.0008730764201389199,0.0008730736531891274,0.0057675121047131486,20.19941734197714,885
added results/inference_v2/run_meta.json +10 −0
@@ -0,0 +1,10 @@
1 +{
2 + "steps": [
3 + "assert",
4 + "cluster",
5 + "lot",
6 + "paired"
7 + ],
8 + "B": 9999,
9 + "elapsed_s": 92
10 +}
\ No newline at end of file
added results/inference_v2/significance_counts_v2.csv +4 −0
@@ -0,0 +1,4 @@
1 +se,n_slopes,n_abs_t_gt_196,critical_t,n_abs_t_gt_crit,n_holm_5pct,n_bh_5pct,median_G_eff_css,min_G_eff_css
2 +hc3,63,58,1.96,58,54,58,,
3 +cluster_zip3,63,40,1.9626481243816396,40,23,36,163.98534085042598,22.031024576681364
4 +cluster_state,63,33,2.0085591121007607,33,17,26,14.229812245284245,3.663564544662053
added results/inference_v2/zip3_missing_cluster.csv +2 −0
@@ -0,0 +1,2 @@
1 +n_listings_missing_zip,n_zip3_prefixes,n_clusters_with_000,max_abs_rel_se_change_after_reassignment,reassigned_to
2 +16,885,886,0.00026019561697609994,119 146 265 283 297 312 367 379 454 545 648 810 864 890 939 982
added results/inference_v2/zip3_missing_reassignment.csv +65 −0
@@ -0,0 +1,65 @@
1 +variable,se_cluster_zip3_with_000,se_cluster_zip3_reassigned,rel_change
2 +const,0.11084286267316142,0.11084313226914856,2.4322358755934914e-06
3 +ln_sqft,0.0161792092342057,0.01617937113644096,1.0006807682483299e-05
4 +bedrooms,0.012662345532639974,0.012662939565514416,4.6913336309684084e-05
5 +bathrooms,0.015687581211757787,0.01568973796250829,0.00013748140783409468
6 +age,0.06504842474426516,0.06504979854654606,2.1119685623460782e-05
7 +age_sq,0.008989950579174545,0.008991441924089044,0.00016589022390767205
8 +stories,0.005807423758672755,0.005807454597652511,5.310268552438657e-06
9 +bath_per_bed,0.011785299316252374,0.011785580240748921,2.3836857173442638e-05
10 +sqft_per_bed,0.0063335917623245554,0.006333474781558496,-1.8469893616313016e-05
11 +ln_lot,0.007551897953306306,0.007551997197614064,1.3141637820179852e-05
12 +has_pool,0.02793098670347089,0.027929725170333493,-4.5166078477265614e-05
13 +has_spa,0.01698996552718787,0.01699006429169212,5.813107983865251e-06
14 +has_basement,0.0214857335273043,0.02148383178415964,-8.851190219993921e-05
15 +has_fireplace,0.011531290143509391,0.011533657515468872,0.00020529983462536094
16 +has_garage,0.010336851345120251,0.01033879387376652,0.00018792266439882077
17 +on_waterfront,0.24417198192494594,0.2441736872145198,6.983969087714925e-06
18 +has_central_air,0.01394817717210821,0.01394982274923529,0.00011797793408963031
19 +has_forced_air,0.01516491952766096,0.01516405484259848,-5.701877025476865e-05
20 +has_hardwood,0.00944439627051184,0.009446853661026412,0.00026019561697609994
21 +parking_spaces,0.004314311335556711,0.004314708156855645,9.197790054304988e-05
22 +luxury_score,0.008549381768231569,0.00854967963167534,3.484034891010346e-05
23 +walk_score,0.012397583396033071,0.01239935293243422,0.000142732365221665
24 +bike_score,0.012152014141512105,0.012152905193339307,7.332544357052662e-05
25 +transit_score,0.01809145703790151,0.018091574987403755,6.519624262235624e-06
26 +avg_school_rating,0.013389621998098181,0.013389743549005077,9.077993905481918e-06
27 +school_count,0.002510806710253446,0.002510804564089567,-8.54770648150982e-07
28 +nearest_school_distance,0.005430923137351372,0.0054309741150743674,9.386566833358145e-06
29 +property_tax_rate,0.014635447354455232,0.014635872637121834,2.9058398851811162e-05
30 +is_condo,0.040180853579487565,0.040180936374243845,2.060552449778541e-06
31 +has_hoa,0.02246088032821604,0.022460874357845953,-2.6581193612873477e-07
32 +ln_hoa,0.012635782200760275,0.012635770526180458,-9.239301240748787e-07
33 +tag_new_construction,0.02278278972506957,0.022781726134127253,-4.6683964306071246e-05
34 +tag_foreclosure,0.022168248586423734,0.022168423623329075,7.895838259885934e-06
35 +sqft_x_age,0.06464022434229162,0.0646427127466216,3.849622050822177e-05
36 +pool_x_south,0.015360411081263122,0.015360250099556287,-1.0480299386661507e-05
37 +waterfront_x_sqft,0.06815329712442604,0.06815367102343607,5.486147050870471e-06
38 +basement_x_north,0.012133956179348968,0.012133597918214216,-2.952550095425277e-05
39 +condo_x_walkscore,0.01558509183297419,0.015585031809819244,-3.851318656988845e-06
40 +age_x_luxury,0.00827906339835743,0.008279311189902728,2.9929900687397648e-05
41 +roof_cat_Metal,0.04349388062002295,0.04349435273069271,1.0854645826618992e-05
42 +roof_cat_Other,0.03934087130984624,0.0393410411421759,4.316943779114624e-06
43 +roof_cat_Shingle_Comp,0.04194027562623559,0.04194061701925363,8.139980316013151e-06
44 +roof_cat_Slate,0.05103160297577027,0.05103228194244627,1.3304827526594565e-05
45 +roof_cat_Tile,0.04466982958354879,0.04467110875965762,2.8636243315949272e-05
46 +roof_cat_Unknown,0.04770444842551859,0.04770498920582322,1.1336056122246418e-05
47 +construction_cat_Concrete_Block,0.03682500870417783,0.03682519158580489,4.966234455672236e-06
48 +construction_cat_Other,0.01291572716024121,0.012915105964390272,-4.8096080323634105e-05
49 +construction_cat_Steel,0.038287788169592996,0.03828786525150359,2.013224432273475e-06
50 +construction_cat_Stone,0.014723890394726274,0.014723897586874764,4.884679454253416e-07
51 +construction_cat_Stucco,0.02606935999504513,0.026069667220825873,1.1784937597347778e-05
52 +construction_cat_Unknown,0.026759383209593708,0.0267618224718596,9.115539946447093e-05
53 +construction_cat_Vinyl,0.011434497065796292,0.011434531715468528,3.0302751434785336e-06
54 +construction_cat_Wood_Frame,0.015193432991140286,0.015194201363075442,5.057263461161199e-05
55 +foundation_cat_Concrete,0.07990891509398913,0.07990883737112651,-9.7264319665058e-07
56 +foundation_cat_Crawlspace,0.08189368995768126,0.08189377501895712,1.038679242615359e-06
57 +foundation_cat_Other,0.08176791319712322,0.08176791168763399,-1.846065489363724e-08
58 +foundation_cat_Pier_Raised,0.08650916880832156,0.0865092095489394,4.7093988309043766e-07
59 +foundation_cat_Slab,0.07981374203303475,0.07981383735744013,1.1943357491883688e-06
60 +foundation_cat_Stone,0.08938436025355444,0.08938524549264738,9.903735848526551e-06
61 +foundation_cat_Unknown,0.07895862665340403,0.07895886488055369,3.017113642478364e-06
62 +region_Northeast,0.05357587517963916,0.053583633940610546,0.00014481818440437166
63 +region_South,0.03385698865833608,0.03385686878415574,-3.540603730312064e-06
64 +region_West,0.047435232182475966,0.04743809230992967,6.029542435270763e-05
65 +lot_missing,0.026876878999291056,0.026877300882287323,1.569687448754742e-05
added scripts/06_inference_v2.py +508 −0
@@ -0,0 +1,508 @@
1 +#!/usr/bin/env python3
2 +# Author: Simon-Pierre Boucher — contact@spboucher.ai
3 +#
4 +"""Inference on the corrected design (November 2026 review, batch L2; items
5 +WP03-DATA-06, EDIT-04, ADVE-05, ADVE-08, ECON-02, THEO-06, ECON-04, THEO-05,
6 +EDIT-07, ADVE-14, ADVE-12, ECON-10, THEO-08, ECON-07, EDIT-09, ECON-08, ADVE-23,
7 +EDIT-02). Refactor of 06_revision_inference.py on the v2 design (trimmed sample).
8 +
9 +Steps (--steps to select; default all)
10 + assert zpid / split-index assertions (DATA-06)
11 + cluster OLS with HC3 / state-CR1 / ZIP3-CR1 SEs; counts over the 62 slopes with
12 + |t| > 1.96, t(G-1), Holm, BH; CSS effective number of clusters; '000' cluster
13 + reported and reassigned to the nearest ZIP3 centroid (SEs compared);
14 + clustered Kennedy effects for the dummies
15 + wild wild cluster bootstrap by state (wildboottest, Rademacher and Webb,
16 + B = 9,999, seed 42) for every slope
17 + states state-by-state OLS with ZIP3-cluster SEs, reversal counts (|t| > 2; state CI
18 + excludes the pooled value), cluster-robust homogeneity test, EB shrinkage
19 + lot lot-size codings on the v2 sample (legacy vs v2, trimmed and untrimmed)
20 + moran Moran's I on the full random test set (cKDTree k = 8, 999 permutations) and
21 + by distance band (0-1, 1-5, 5-25 km) for OLS / LightGBM / XGBoost residuals
22 + (stored v1.2 predictions and, when present, v2 predictions from batch N2)
23 + conley Conley HAC (Bartlett, 25 / 50 / 100 km) on a 100,000-listing subsample
24 + paired ZIP3-cluster SE of paired squared-error differences on the random split
25 +
26 +Outputs: results/inference_v2/*.csv
27 +
28 +Usage: python scripts/06_inference_v2.py [--steps assert cluster wild states lot moran conley paired]
29 + [--B 9999] [--workers 8]
30 +"""
31 +
32 +import argparse
33 +import json
34 +import sys
35 +import time
36 +import warnings
37 +from pathlib import Path
38 +
39 +import numpy as np
40 +import pandas as pd
41 +import statsmodels.api as sm
42 +from scipy import sparse, stats
43 +from scipy.spatial import cKDTree
44 +from sklearn.metrics import r2_score
45 +from sklearn.model_selection import train_test_split
46 +
47 +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
48 +from wp3 import config, data, design # noqa: E402
49 +
50 +OUT = config.RESULTS_DIR / "inference_v2"
51 +KEY24 = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq", "ln_lot", "lot_missing", "has_pool",
52 + "has_garage", "on_waterfront", "has_central_air", "has_hardwood", "luxury_score",
53 + "walk_score", "bike_score", "transit_score", "avg_school_rating", "property_tax_rate",
54 + "is_condo", "has_hoa", "tag_new_construction", "tag_foreclosure",
55 + "region_Northeast", "region_South", "region_West"]
56 +REPORT13 = ["ln_sqft", "bedrooms", "bathrooms", "age", "ln_lot", "has_pool", "has_garage",
57 + "has_central_air", "has_hardwood", "luxury_score", "tag_foreclosure", "is_condo", "has_hoa"]
58 +DUMMIES = ["on_waterfront", "has_pool", "has_garage", "has_central_air", "has_hardwood", "is_condo",
59 + "tag_foreclosure", "tag_new_construction", "has_hoa", "has_basement", "has_spa", "has_fireplace", "lot_missing"]
60 +
61 +
62 +def holm(p):
63 + p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 0.0
64 + for i, j in enumerate(order):
65 + v = max(min(1.0, (m - i) * p[j]), prev); adj[j] = v; prev = v
66 + return adj
67 +
68 +
69 +def bh(p):
70 + p = np.asarray(p, float); m = len(p); order = np.argsort(p); adj = np.empty(m); prev = 1.0
71 + for i in range(m - 1, -1, -1):
72 + j = order[i]; v = min(min(1.0, m / (i + 1) * p[j]), prev); adj[j] = v; prev = v
73 + return adj
74 +
75 +
76 +def to_xy_km(lat, lon):
77 + """Equirectangular projection (km), adequate for neighbour searches below ~100 km."""
78 + lat0 = np.deg2rad(np.mean(lat))
79 + return np.column_stack([np.deg2rad(lon) * 6371.0 * np.cos(lat0), np.deg2rad(lat) * 6371.0])
80 +
81 +
82 +# ----------------------------------------------------------------------------- assert
83 +def step_assert():
84 + sample = pd.read_parquet(config.DATA_PROCESSED / "analytical_sample_v2.parquet", columns=["zpid"])
85 + d = pd.read_parquet(design.DESIGN_V2, columns=["zpid"])
86 + assert (sample["zpid"].values == d["zpid"].values).all(), "design_v2 / analytical_sample_v2 zpid order differs"
87 + rows = {"design_vs_sample_zpid_order": True}
88 + if config.ANALYTICAL_SAMPLE.exists() and config.MODEL_DATA.exists():
89 + pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)
90 + md = data.load_model_data()
91 + assert (md["df_clean"]["zpid"].values == pk["zpid"].values).all(), "model_data.pkl / analytical_sample.pkl zpid order differs"
92 + assert (pk["zpid"].values == d["zpid"].values).all()
93 + rows["model_data_vs_pickle_zpid_order"] = True
94 + try:
95 + ml = data.load_result("ml_results.pkl")
96 + idx_tr, idx_te = train_test_split(np.arange(len(pk)), test_size=0.2, random_state=config.RANDOM_STATE)
97 + assert (ml["X_test"].index.values == idx_te).all(), "ml_results X_test index != train_test_split(seed 42)"
98 + assert ml["X_test"].index.isin(pk.index).all()
99 + rows["ml_results_test_index_equals_seed42_split"] = True
100 + except FileNotFoundError:
101 + rows["ml_results_test_index_equals_seed42_split"] = "ml_results.pkl absent"
102 + pd.Series(rows).to_csv(OUT / "assertions.csv", header=["passed"])
103 + print("assertions:", rows)
104 +
105 +
106 +# ----------------------------------------------------------------------------- cluster
107 +def css_effective_clusters(X, groups):
108 + """Carter-Schnepel-Steigerwald (2017) feasible effective number of clusters per coefficient
109 + (homoskedastic within-cluster errors, Omega_g = I): gamma_g = e_k' A X_g'X_g A e_k."""
110 + A = np.linalg.inv(X.T @ X)
111 + codes, uniq = pd.factorize(groups)
112 + G = len(uniq); K = X.shape[1]
113 + gam = np.zeros((G, K))
114 + for g in range(G):
115 + Xg = X[codes == g]
116 + M = A @ (Xg.T @ Xg) @ A
117 + gam[g] = np.diag(M)
118 + gbar = gam.mean(axis=0)
119 + Gamma = ((gam - gbar) ** 2 / gbar ** 2).mean(axis=0)
120 + return G / (1 + Gamma)
121 +
122 +
123 +def step_cluster():
124 + X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)
125 + Xc = sm.add_constant(X, has_constant="add")
126 + names = list(Xc.columns)
127 + st = pd.factorize(info["address_state"].values)[0]
128 + z3 = pd.factorize(info["zip3"].values)[0]
129 + base = sm.OLS(y, Xc.values)
130 + r_hc3 = base.fit(cov_type="HC3")
131 + r_st = base.fit(cov_type="cluster", cov_kwds={"groups": st})
132 + r_z3 = base.fit(cov_type="cluster", cov_kwds={"groups": z3})
133 + out = pd.DataFrame({"coef": r_hc3.params, "se_hc3": r_hc3.bse, "se_cluster_state": r_st.bse,
134 + "se_cluster_zip3": r_z3.bse}, index=names)
135 + out["ratio_zip3_hc3"] = out.se_cluster_zip3 / out.se_hc3; out["ratio_state_hc3"] = out.se_cluster_state / out.se_hc3
136 + for k in ["hc3", "cluster_state", "cluster_zip3"]:
137 + out[f"t_{k}"] = out.coef / out[f"se_{k}"]
138 + out["p_state_t50"] = 2 * (1 - stats.t.cdf(out.t_cluster_state.abs(), 50))
139 + out["p_zip3_t885"] = 2 * (1 - stats.t.cdf(out.t_cluster_zip3.abs(), len(np.unique(z3)) - 1))
140 + out["p_hc3_normal"] = 2 * (1 - stats.norm.cdf(out.t_hc3.abs()))
141 + slopes = out.index != "const"
142 + for k, pcol in [("state", "p_state_t50"), ("zip3", "p_zip3_t885"), ("hc3", "p_hc3_normal")]:
143 + out.loc[slopes, f"p_holm_{k}"] = holm(out.loc[slopes, pcol]); out.loc[slopes, f"p_bh_{k}"] = bh(out.loc[slopes, pcol])
144 + # CSS effective clusters (state)
145 + out["G_eff_css_state"] = css_effective_clusters(Xc.values, info["address_state"].values)
146 + out["G_eff_css_zip3"] = css_effective_clusters(Xc.values, info["zip3"].values)
147 + out["n_clusters_state"] = len(np.unique(st)); out["n_clusters_zip3"] = len(np.unique(z3))
148 + out.index.name = "variable"
149 + out.to_csv(OUT / "cluster_se_v2.csv")
150 + s = out[slopes]
151 + counts = []
152 + for k, G in [("hc3", None), ("cluster_zip3", len(np.unique(z3))), ("cluster_state", 51)]:
153 + t = s[f"t_{k}"].abs()
154 + crit = stats.t.ppf(0.975, G - 1) if G else 1.96
155 + pc = {"hc3": "p_hc3_normal", "cluster_zip3": "p_zip3_t885", "cluster_state": "p_state_t50"}[k]
156 + kk = k.replace("cluster_", "")
157 + counts.append({"se": k, "n_slopes": int(slopes.sum()), "n_abs_t_gt_196": int((t > 1.96).sum()),
158 + "critical_t": crit, "n_abs_t_gt_crit": int((t > crit).sum()),
159 + "n_holm_5pct": int((s[f"p_holm_{kk}"] < 0.05).sum()), "n_bh_5pct": int((s[f"p_bh_{kk}"] < 0.05).sum()),
160 + "median_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].median()) if G else np.nan,
161 + "min_G_eff_css": float(s["G_eff_css_state" if kk == "state" else "G_eff_css_zip3"].min()) if G else np.nan})
162 + pd.DataFrame(counts).to_csv(OUT / "significance_counts_v2.csv", index=False)
163 + print(pd.DataFrame(counts).to_string(index=False))
164 + # '000' cluster: size and nearest-ZIP3 reassignment
165 + miss = info["zip3"].values == "000"
166 + cent = info.loc[~miss].groupby("zip3")[["latitude", "longitude"]].mean()
167 + tree = cKDTree(to_xy_km(cent.latitude.values, cent.longitude.values))
168 + _, j = tree.query(to_xy_km(info.loc[miss, "latitude"].values, info.loc[miss, "longitude"].values))
169 + z3b = info["zip3"].values.copy(); z3b[miss] = cent.index.values[j]
170 + r_z3b = base.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3b)[0]})
171 + cmp = pd.DataFrame({"se_cluster_zip3_with_000": r_z3.bse, "se_cluster_zip3_reassigned": r_z3b.bse}, index=names)
172 + cmp["rel_change"] = cmp.se_cluster_zip3_reassigned / cmp.se_cluster_zip3_with_000 - 1
173 + cmp.index.name = "variable"; cmp.to_csv(OUT / "zip3_missing_reassignment.csv")
174 + pd.DataFrame([{"n_listings_missing_zip": int(miss.sum()), "n_zip3_prefixes": int(len(cent)),
175 + "n_clusters_with_000": int(len(np.unique(z3))), "max_abs_rel_se_change_after_reassignment": float(cmp.rel_change.abs().max()),
176 + "reassigned_to": " ".join(sorted(set(cent.index.values[j])))}]).to_csv(OUT / "zip3_missing_cluster.csv", index=False)
177 + # Kennedy with clustered variance (natural units)
178 + Xu, yu, infou, _ = design.load_design("v2", trimmed=True, standardize=False)
179 + Xcu = sm.add_constant(Xu, has_constant="add")
180 + bu = sm.OLS(yu, Xcu.values)
181 + ru = {"hc3": bu.fit(cov_type="HC3"), "cluster_zip3": bu.fit(cov_type="cluster", cov_kwds={"groups": z3}),
182 + "cluster_state": bu.fit(cov_type="cluster", cov_kwds={"groups": st})}
183 + kr = []
184 + for v in DUMMIES:
185 + jv = list(Xcu.columns).index(v)
186 + for k, r in ru.items():
187 + b, se_ = r.params[jv], r.bse[jv]
188 + kr.append({"variable": v, "se": k, "coef": b, "se_value": se_, "naive_pct": 100 * (np.exp(b) - 1),
189 + "kennedy_pct": 100 * (np.exp(b - 0.5 * se_ ** 2) - 1), "t": b / se_})
190 + pd.DataFrame(kr).to_csv(OUT / "kennedy_clustered_v2.csv", index=False)
191 + unstd = pd.DataFrame({"coef": ru["hc3"].params, "se_hc3": ru["hc3"].bse, "se_cluster_zip3": ru["cluster_zip3"].bse,
192 + "se_cluster_state": ru["cluster_state"].bse}, index=list(Xcu.columns))
193 + unstd.index.name = "variable"; unstd.to_csv(OUT / "ols_v2_unstandardized_cluster_se.csv")
194 +
195 +
196 +# ----------------------------------------------------------------------------- wild
197 +def step_wild(B, workers):
198 + from wildboottest.wildboottest import wildboottest
199 + X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)
200 + Xc = sm.add_constant(X, has_constant="add")
201 + model = sm.OLS(y, Xc)
202 + cl = pd.Series(pd.factorize(info["address_state"].values)[0])
203 + rows = []
204 + t0 = time.time()
205 + for wt in ["rademacher", "webb"]:
206 + for v in [c for c in Xc.columns if c != "const"]:
207 + res = wildboottest(model, cluster=cl, B=B, param=v, weights_type=wt, seed=config.RANDOM_STATE, show=False)
208 + r = res.iloc[0].to_dict()
209 + rows.append({"variable": v, "weights": wt, **{k: r[k] for k in r if k in ("statistic", "p-value", "param")}, "B": B})
210 + print(f" wild bootstrap ({wt}) done in {time.time() - t0:.0f}s", flush=True)
211 + out = pd.DataFrame(rows).rename(columns={"p-value": "p_wild", "statistic": "t_cr1"})
212 + wide = out.pivot(index="variable", columns="weights", values="p_wild").rename(columns=lambda c: f"p_wild_{c}")
213 + tt = out[out.weights == "webb"].set_index("variable")["t_cr1"]
214 + wide["t_cr1"] = tt
215 + wide["p_holm_webb"] = holm(wide["p_wild_webb"]); wide["p_bh_webb"] = bh(wide["p_wild_webb"])
216 + wide["p_holm_rademacher"] = holm(wide["p_wild_rademacher"]); wide["p_bh_rademacher"] = bh(wide["p_wild_rademacher"])
217 + wide.to_csv(OUT / "wild_bootstrap_state.csv")
218 + cnt = {"n_slopes": len(wide), "n_sig_5pct_webb": int((wide.p_wild_webb < 0.05).sum()),
219 + "n_sig_5pct_rademacher": int((wide.p_wild_rademacher < 0.05).sum()),
220 + "n_sig_5pct_webb_holm": int((wide.p_holm_webb < 0.05).sum()), "n_sig_5pct_webb_bh": int((wide.p_bh_webb < 0.05).sum()),
221 + "n_sig_5pct_rademacher_holm": int((wide.p_holm_rademacher < 0.05).sum()), "B": B}
222 + pd.DataFrame([cnt]).to_csv(OUT / "wild_bootstrap_counts.csv", index=False)
223 + print(cnt)
224 +
225 +
226 +# ----------------------------------------------------------------------------- states
227 +def step_states():
228 + Xu, y, info, _ = design.load_design("v2", trimmed=True, standardize=False)
229 + cols = [c for c in Xu.columns if not c.startswith("region_")]
230 + states = info["address_state"].values; z3 = info["zip3"].values
231 + idx_tr, idx_te = train_test_split(np.arange(788_842), test_size=0.2, random_state=config.RANDOM_STATE)
232 + # map the full-sample split onto the trimmed sample via zpid
233 + full = pd.read_parquet(design.DESIGN_V2, columns=["zpid"])
234 + test_z = set(full["zpid"].values[idx_te])
235 + is_test = info["zpid"].isin(test_z).values
236 + rows, se_rows, dropped = [], [], []
237 + pred_inter = np.full(len(y), np.nan)
238 + for s_ in sorted(np.unique(states)):
239 + m = states == s_
240 + Xs = Xu.loc[m, cols]
241 + keep = [c for c in cols if Xs[c].nunique() > 1]
242 + dropped.append({"state": s_, "n": int(m.sum()), "dropped_constant_columns": " ".join(c for c in cols if c not in keep)})
243 + Xc = sm.add_constant(Xs[keep].values, has_constant="add")
244 + g = pd.factorize(z3[m])[0]
245 + r = sm.OLS(y[m], Xc).fit(cov_type="cluster", cov_kwds={"groups": g}) if len(np.unique(g)) > 1 else sm.OLS(y[m], Xc).fit(cov_type="HC3")
246 + cond = np.linalg.cond(Xc)
247 + row = {"state": s_, "n": int(m.sum()), "n_zip3": int(len(np.unique(g))), "r2": r.rsquared, "cond_number": cond,
248 + "median_price": float(np.exp(np.median(y[m])))}
249 + for v in REPORT13:
250 + if v in keep:
251 + j = keep.index(v) + 1
252 + row[v] = r.params[j]; row[f"se_{v}"] = r.bse[j]
253 + else:
254 + row[v] = np.nan; row[f"se_{v}"] = np.nan
255 + rows.append(row)
256 + tr = m & ~is_test; te = m & is_test
257 + if tr.sum() > len(keep) + 5 and te.sum() > 0:
258 + from sklearn.linear_model import LinearRegression
259 + pred_inter[te] = LinearRegression().fit(Xu.loc[tr, keep].values, y[tr]).predict(Xu.loc[te, keep].values)
260 + het = pd.DataFrame(rows); het.to_csv(OUT / "state_heterogeneity_v2.csv", index=False)
261 + pd.DataFrame(dropped).to_csv(OUT / "state_regressions_dropped_columns.csv", index=False)
262 + # pooled state-FE reference with ZIP3-cluster SEs
263 + sfe = pd.get_dummies(pd.Series(states), drop_first=True).astype(float).values
264 + Xp = np.column_stack([np.ones(len(y)), Xu[cols].values, sfe])
265 + pooled = sm.OLS(y, Xp).fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(z3)[0]})
266 + summ = []
267 + for v in REPORT13:
268 + jp = cols.index(v) + 1
269 + b0, s0 = pooled.params[jp], pooled.bse[jp]
270 + b = het[v].values; se_ = het[f"se_{v}"].values; ok = ~np.isnan(b) & ~np.isnan(se_) & (se_ > 0)
271 + b, se_ = b[ok], se_[ok]
272 + w = 1 / se_ ** 2
273 + b_pw = np.sum(w * b) / np.sum(w)
274 + # cluster-robust homogeneity (states independent; ZIP3 clusters nested in states)
275 + Q = float(np.sum(w * (b - b_pw) ** 2)); dfq = len(b) - 1
276 + tau2 = max(0.0, (Q - dfq) / (np.sum(w) - np.sum(w ** 2) / np.sum(w))) # DerSimonian-Laird
277 + summ.append({"variable": v, "pooled_state_fe": b0, "pooled_se_zip3": s0, "n_states": int(len(b)),
278 + "precision_weighted_mean": b_pw, "listing_weighted_mean": np.average(het[v].values[ok], weights=het["n"].values[ok]),
279 + "p10": np.percentile(b, 10), "median": np.median(b), "p90": np.percentile(b, 90), "min": b.min(), "max": b.max(),
280 + "n_sign_differs_from_pooled": int((np.sign(b) != np.sign(b0)).sum()),
281 + "n_sign_differs_and_abs_t_gt_2": int(((np.sign(b) != np.sign(b0)) & (np.abs(b / se_) > 2)).sum()),
282 + "n_state_ci_excludes_pooled": int(((b - 1.96 * se_ > b0) | (b + 1.96 * se_ < b0)).sum()),
283 + "n_abs_t_gt_2": int((np.abs(b / se_) > 2).sum()),
284 + "homogeneity_Q": Q, "homogeneity_df": dfq, "homogeneity_p": float(1 - stats.chi2.cdf(Q, dfq)),
285 + "tau2_dl": tau2, "tau_dl": np.sqrt(tau2), "I2_pct": max(0.0, 100 * (Q - dfq) / Q) if Q > 0 else 0.0})
286 + summ = pd.DataFrame(summ); summ.to_csv(OUT / "state_heterogeneity_summary_v2.csv", index=False)
287 + print(summ[["variable", "pooled_state_fe", "p10", "median", "p90", "n_sign_differs_from_pooled", "n_sign_differs_and_abs_t_gt_2",
288 + "n_state_ci_excludes_pooled", "homogeneity_p", "tau_dl"]].round(3).to_string(index=False))
289 + # empirical-Bayes shrinkage for Figure 3 (normal-normal, tau2 DL)
290 + eb = []
291 + for v in REPORT13:
292 + t2 = float(summ.set_index("variable").loc[v, "tau2_dl"]); mu = float(summ.set_index("variable").loc[v, "precision_weighted_mean"])
293 + for _, r in het.iterrows():
294 + b, s_ = r[v], r[f"se_{v}"]
295 + if np.isnan(b) or np.isnan(s_) or s_ <= 0:
296 + continue
297 + shrink = t2 / (t2 + s_ ** 2)
298 + post = shrink * b + (1 - shrink) * mu; pv = shrink * s_ ** 2
299 + eb.append({"variable": v, "state": r.state, "n": r.n, "estimate": b, "se_zip3": s_, "ci_low": b - 1.96 * s_, "ci_high": b + 1.96 * s_,
300 + "eb_posterior_mean": post, "eb_posterior_sd": np.sqrt(pv), "eb_low": post - 1.96 * np.sqrt(pv), "eb_high": post + 1.96 * np.sqrt(pv),
301 + "shrinkage_factor": shrink})
302 + pd.DataFrame(eb).to_csv(OUT / "state_heterogeneity_eb.csv", index=False)
303 + # pooled vs state FE vs fully interacted on the random test split (trimmed)
304 + from sklearn.linear_model import LinearRegression
305 + te = is_test; tr = ~is_test
306 + p_pooled = LinearRegression().fit(Xu.loc[tr, cols].values, y[tr]).predict(Xu.loc[te, cols].values)
307 + gtr = pd.Series(states[tr]); Xd = Xu.loc[tr, cols].values - pd.DataFrame(Xu.loc[tr, cols].values).groupby(gtr.values).transform("mean").values
308 + yd = y[tr] - pd.Series(y[tr]).groupby(gtr.values).transform("mean").values
309 + beta, *_ = np.linalg.lstsq(Xd, yd, rcond=None)
310 + alpha = pd.Series(y[tr] - Xu.loc[tr, cols].values @ beta).groupby(gtr.values).mean()
311 + p_fe = Xu.loc[te, cols].values @ beta + pd.Series(states[te]).map(alpha).fillna(alpha.mean()).values
312 + ok = ~np.isnan(pred_inter[te])
313 + comp = pd.DataFrame([{"model": "Pooled OLS, no geographic controls", "r2_test": r2_score(y[te], p_pooled)},
314 + {"model": "Pooled OLS + 50 state fixed effects", "r2_test": r2_score(y[te], p_fe)},
315 + {"model": "State-by-state OLS (fully interacted)", "r2_test": r2_score(y[te][ok], pred_inter[te][ok])}])
316 + comp["n_train"] = int(tr.sum()); comp["n_test"] = int(te.sum())
317 + comp.to_csv(OUT / "pooled_vs_interacted_v2.csv", index=False); print(comp.to_string(index=False))
318 +
319 +
320 +# ----------------------------------------------------------------------------- lot
321 +def step_lot():
322 + rows = []
323 + for trimmed in [False, True]:
324 + Xl, y, info, _ = design.load_design("legacy", trimmed=trimmed, standardize=False, lot_indicator=False)
325 + Xv, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=False)
326 + Xvi, _, _, _ = design.load_design("v2", trimmed=trimmed, standardize=False, lot_indicator=True)
327 + lm = info["lot_missing"].values == 1
328 + def fit(X, mask, label):
329 + Xc = sm.add_constant(X.values[mask], has_constant="add")
330 + r = sm.OLS(y[mask], Xc).fit(cov_type="HC3")
331 + names = list(X.columns)
332 + return {"sample": "trimmed (783,872)" if trimmed else "all (788,842)", "specification": label, "N": int(mask.sum()),
333 + "R2": r.rsquared, "beta_ln_lot": r.params[names.index("ln_lot") + 1], "se_ln_lot": r.bse[names.index("ln_lot") + 1],
334 + "beta_ln_sqft": r.params[names.index("ln_sqft") + 1],
335 + "beta_lot_missing": r.params[names.index("lot_missing") + 1] if "lot_missing" in names else np.nan}
336 + allm = np.ones(len(y), bool)
337 + rows.append(fit(Xl, allm, "v1.1 coding: ln(lot) = 0 when missing, zero, 1 sq ft or negative"))
338 + rows.append(fit(Xv, allm, "v2: ln(state-median lot) when missing/zero/<= 1"))
339 + rows.append(fit(Xvi, allm, "v2 + lot-missing indicator (the design used throughout)"))
340 + rows.append(fit(Xl, ~lm, "Complete cases (raw lot > 1 sq ft)"))
341 + p995 = info["lot_size_sqft"].quantile(0.995)
342 + Xw = Xvi.copy(); lot_v2 = np.exp(Xvi["ln_lot"].values); Xw["ln_lot"] = np.log(np.minimum(lot_v2, p995))
343 + rows.append(fit(Xw, allm, "v2 + indicator, lot winsorised at p99.5"))
344 + out = pd.DataFrame(rows); out.to_csv(OUT / "lot_recoding_v2.csv", index=False); print(out.round(4).to_string(index=False))
345 +
346 +
347 +# ----------------------------------------------------------------------------- moran
348 +def moran_knn(e, xy, k=8, permutations=999, seed=42):
349 + from esda.moran import Moran
350 + from libpysal.weights import WSP
351 + tree = cKDTree(xy)
352 + d, j = tree.query(xy, k=k + 1)
353 + n = len(xy)
354 + rows_ = np.repeat(np.arange(n), k); cols_ = j[:, 1:].ravel()
355 + W = sparse.csr_matrix((np.ones(n * k) / k, (rows_, cols_)), shape=(n, n))
356 + w = WSP(W).to_W()
357 + np.random.seed(seed)
358 + mi = Moran(e, w, transformation="r", permutations=permutations)
359 + return {"I": mi.I, "E_I": mi.EI, "z_norm": mi.z_norm, "p_norm": mi.p_norm, "z_sim": mi.z_sim, "p_sim": mi.p_sim,
360 + "median_kth_neighbour_km": float(np.median(d[:, k]))}
361 +
362 +
363 +def moran_band(e, pairs_i, pairs_j, n):
364 + """Moran's I for a symmetric binary band weight matrix given as pair lists; normal-approx z (randomisation)."""
365 + z = e - e.mean()
366 + S0 = 2 * len(pairs_i)
367 + num = 2 * np.sum(z[pairs_i] * z[pairs_j])
368 + I = (n / S0) * num / np.sum(z ** 2)
369 + EI = -1 / (n - 1)
370 + deg = np.bincount(np.concatenate([pairs_i, pairs_j]), minlength=n).astype(float)
371 + S1 = 0.5 * np.sum((2.0) ** 2) * len(pairs_i) # each symmetric pair contributes (w_ij + w_ji)^2 = 4, times 1/2
372 + S2 = np.sum((2 * deg) ** 2)
373 + b2 = n * np.sum(z ** 4) / np.sum(z ** 2) ** 2
374 + VI = (n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2) - b2 * ((n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)) / \
375 + ((n - 1) * (n - 2) * (n - 3) * S0 ** 2) - EI ** 2
376 + zI = (I - EI) / np.sqrt(VI)
377 + return {"I": I, "E_I": EI, "z_rand": zI, "p_rand": float(2 * (1 - stats.norm.cdf(abs(zI)))), "n_pairs": int(len(pairs_i)),
378 + "share_units_with_neighbour": float((deg > 0).mean())}
379 +
380 +
381 +def step_moran(permutations):
382 + rows = []
383 + sets = []
384 + # v1.2 stored predictions (random test set, legacy design)
385 + try:
386 + ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)
387 + idx = ml["X_test"].index.values
388 + xy = to_xy_km(pk.loc[idx, "latitude"].values, pk.loc[idx, "longitude"].values)
389 + y = ml["y_test"].values
390 + for name, p in [("OLS", ml["y_pred_ols"]), ("LightGBM", ml["y_pred_lgb"]), ("XGBoost", ml["y_pred_xgb"])]:
391 + sets.append(("v1.2 stored predictions, random test set", name, y - p, xy))
392 + except FileNotFoundError:
393 + pass
394 + # v2 predictions from batch N2 when available
395 + prim = config.RESULTS_DIR / "primary"
396 + d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "latitude", "longitude"]).set_index("zpid")
397 + for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"),
398 + ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"),
399 + ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]:
400 + f = prim / pat
401 + if f.exists():
402 + p = pd.read_parquet(f)
403 + xy = to_xy_km(d.loc[p.zpid, "latitude"].values, d.loc[p.zpid, "longitude"].values)
404 + sets.append(("v2 predictions (N2), random test set, trimmed", name, p["y"].values - p["pred"].values, xy))
405 + bands = [(0, 1), (1, 5), (5, 25)]
406 + for label, name, e, xy in sets:
407 + t0 = time.time()
408 + r = moran_knn(e, xy, k=8, permutations=permutations)
409 + rows.append({"predictions": label, "model": name, "weights": "kNN k=8 row-standardised", "n": len(e), **r})
410 + print(f" {label} {name} kNN I={r['I']:.4f} z_sim={r['z_sim']:.1f} ({time.time() - t0:.0f}s)", flush=True)
411 + tree = cKDTree(xy)
412 + prev = None
413 + for lo, hi in bands:
414 + pairs = tree.query_pairs(hi, output_type="ndarray")
415 + if lo > 0:
416 + dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)
417 + pairs = pairs[dd > lo]
418 + rb = moran_band(e, pairs[:, 0], pairs[:, 1], len(e))
419 + rows.append({"predictions": label, "model": name, "weights": f"distance band {lo}-{hi} km (binary)", "n": len(e), **rb})
420 + print(f" band {lo}-{hi} km: I={rb['I']:.4f} z={rb['z_rand']:.1f} pairs={rb['n_pairs']:,}", flush=True)
421 + out = pd.DataFrame(rows); out.to_csv(OUT / "moran_full_test_v2.csv", index=False)
422 +
423 +
424 +# ----------------------------------------------------------------------------- conley
425 +def step_conley(n_sub=100_000, cutoffs=(25, 50, 100)):
426 + X, y, info, ex = design.load_design("v2", trimmed=True, standardize=True)
427 + rng = np.random.RandomState(config.RANDOM_STATE)
428 + idx = np.sort(rng.choice(len(y), n_sub, replace=False))
429 + Xc = sm.add_constant(X.iloc[idx], has_constant="add").values; yy = y[idx]
430 + names = ["const"] + list(X.columns)
431 + r = sm.OLS(yy, Xc).fit(cov_type="HC0")
432 + e = r.resid
433 + A = np.linalg.inv(Xc.T @ Xc)
434 + xy = to_xy_km(info["latitude"].values[idx], info["longitude"].values[idx])
435 + tree = cKDTree(xy)
436 + out = pd.DataFrame({"coef": r.params, "se_hc0": r.bse}, index=names)
437 + st = pd.factorize(info["address_state"].values[idx])[0]
438 + out["se_cluster_state"] = sm.OLS(yy, Xc).fit(cov_type="cluster", cov_kwds={"groups": st}).bse
439 + Xe = Xc * e[:, None]
440 + for c in cutoffs:
441 + t0 = time.time()
442 + pairs = tree.query_pairs(c, output_type="ndarray")
443 + dd = np.linalg.norm(xy[pairs[:, 0]] - xy[pairs[:, 1]], axis=1)
444 + k = 1 - dd / c # Bartlett kernel
445 + n = len(yy)
446 + K = sparse.coo_matrix((np.concatenate([k, k]), (np.concatenate([pairs[:, 0], pairs[:, 1]]), np.concatenate([pairs[:, 1], pairs[:, 0]]))), shape=(n, n)).tocsr()
447 + M = Xe.T @ Xe + Xe.T @ (K @ Xe)
448 + V = A @ M @ A
449 + out[f"se_conley_{c}km"] = np.sqrt(np.diag(V))
450 + print(f" Conley {c} km: {len(pairs):,} pairs ({time.time() - t0:.0f}s)", flush=True)
451 + for c in cutoffs:
452 + out[f"ratio_conley_{c}km_hc0"] = out[f"se_conley_{c}km"] / out["se_hc0"]
453 + out["n_subsample"] = n_sub; out.index.name = "variable"
454 + out.to_csv(OUT / "conley_hac_v2.csv")
455 + print(out.loc[[v for v in KEY24 if v in out.index], ["coef", "se_hc0", "se_cluster_state"] + [f"se_conley_{c}km" for c in cutoffs]].round(4).to_string())
456 +
457 +
458 +# ----------------------------------------------------------------------------- paired
459 +def step_paired():
460 + prim = config.RESULTS_DIR / "primary"
461 + files = {p.name: p for p in prim.glob("predictions_random_*.parquet")}
462 + if not files:
463 + print(" no v2 random-split predictions yet (batch N2) — using stored v1.2 predictions")
464 + ml = data.load_result("ml_results.pkl"); pk = pd.read_pickle(config.ANALYTICAL_SAMPLE)
465 + z3 = pk.loc[ml["X_test"].index, "address_zipcode"].fillna("000").astype(str).str[:3].values
466 + y = ml["y_test"].values; preds = {"OLS": ml["y_pred_ols"], "LightGBM": ml["y_pred_lgb"], "XGBoost": ml["y_pred_xgb"]}; label = "v1.2 stored"
467 + else:
468 + d = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "zip3"]).set_index("zpid")
469 + preds = {}
470 + for name, pat in [("OLS", "predictions_random_OLS_62_regressors_+_lot_indicator.parquet"),
471 + ("LightGBM", "predictions_random_LightGBM_1000_rounds_early_stopping_seed42.parquet"),
472 + ("XGBoost", "predictions_random_XGBoost_1000_rounds_early_stopping_seed42.parquet")]:
473 + p = pd.read_parquet(prim / pat); preds[name] = p["pred"].values; y = p["y"].values; zp = p["zpid"].values
474 + z3 = d.loc[zp, "zip3"].values; label = "v2 (N2)"
475 + g = pd.factorize(z3)[0]
476 + rows = []
477 + for a, c in [("XGBoost", "OLS"), ("XGBoost", "LightGBM"), ("LightGBM", "OLS")]:
478 + dl = (y - preds[c]) ** 2 - (y - preds[a]) ** 2
479 + r = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="cluster", cov_kwds={"groups": g})
480 + r0 = sm.OLS(dl, np.ones(len(dl))).fit(cov_type="HC0")
481 + rows.append({"predictions": label, "quantity": f"{c} minus {a} squared-error difference", "mean": float(dl.mean()),
482 + "se_iid": float(dl.std(ddof=1) / np.sqrt(len(dl))), "se_hc0": float(r0.bse[0]), "se_cluster_zip3": float(r.bse[0]),
483 + "t_cluster_zip3": float(r.params[0] / r.bse[0]), "n_clusters": int(len(np.unique(g)))})
484 + pd.DataFrame(rows).to_csv(OUT / "paired_loss_cluster_se.csv", index=False); print(pd.DataFrame(rows).round(5).to_string(index=False))
485 +
486 +
487 +def main():
488 + ap = argparse.ArgumentParser(description=__doc__)
489 + ap.add_argument("--steps", nargs="+", default=["assert", "cluster", "wild", "states", "lot", "moran", "conley", "paired"])
490 + ap.add_argument("--B", type=int, default=9999)
491 + ap.add_argument("--workers", type=int, default=8)
492 + ap.add_argument("--permutations", type=int, default=999)
493 + args = ap.parse_args()
494 + OUT.mkdir(parents=True, exist_ok=True)
495 + warnings.filterwarnings("ignore")
496 + import os
497 + os.environ.setdefault("OMP_NUM_THREADS", str(args.workers))
498 + t0 = time.time()
499 + for s in args.steps:
500 + print(f"== {s} ({time.time() - t0:.0f}s)", flush=True)
501 + {"assert": step_assert, "cluster": step_cluster, "wild": lambda: step_wild(args.B, args.workers), "states": step_states,
502 + "lot": step_lot, "moran": lambda: step_moran(args.permutations), "conley": step_conley, "paired": step_paired}[s]()
503 + (OUT / "run_meta.json").write_text(json.dumps({"steps": args.steps, "B": args.B, "elapsed_s": round(time.time() - t0)}, indent=1))
504 + print(f"Done in {time.time() - t0:.0f}s")
505 +
506 +
507 +if __name__ == "__main__":
508 + main()
added scripts/09_geo_inference.py +299 −0
@@ -0,0 +1,299 @@
1 +#!/usr/bin/env python3
2 +# Author: Simon-Pierre Boucher — contact@spboucher.ai
3 +#
4 +"""Inference on the geographic designs from the out-of-fold predictions of
5 +05_geo_designs.py (November 2026 review, batch L2; items WP03-ADVE-17, EDIT-01,
6 +THEO-03, THEO-01, THEO-02, ECON-05, EDIT-14, THEO-12, ADVE-22, THEO-14, EDIT-13,
7 +THEO-15, ADVE-03, ECON-06, EDIT-06).
8 +
9 +From results/geo/predictions/<design>_<model>_<variant>_seed<k>.parquet:
10 + gap ZIP3-cluster bootstrap (1,000 draws, seed 42) of R2(OLS) - R2(XGB nested),
11 + - R2(XGB 1,000 rounds), - R2(XGB 100 rounds), - R2(LightGBM nested) for every
12 + design (pooled out-of-fold R2); seed dispersion (mean, SD, min, max over seeds)
13 + loso per-state paired differences (51 states), sign test and Wilcoxon; LOSO
14 + summaries (median, listing-weighted mean, p10-p90, states > 0.5)
15 + decomposition SSE within / between (level) shares and oracle within-group R2 per design
16 + feasible feasible level correction: level estimated from k in {5, 20, 50, 100, 200}
17 + random listings of the held-out state, 100 repetitions (LOSO and ten-state)
18 + two_by_two region in/out x rounds {nested, 100, 1,000} under ten-state and LOSO
19 + duplicates random-split R2 without the test rows that have an exact duplicate in train
20 + range the OLS-vs-XGB gap as a range across designs and seeds (cluster A headline)
21 +
22 +Outputs: results/geo_inference/*.csv
23 +
24 +Usage: python scripts/09_geo_inference.py [--draws 1000]
25 +"""
26 +
27 +import argparse
28 +import sys
29 +import time
30 +from pathlib import Path
31 +
32 +import numpy as np
33 +import pandas as pd
34 +from scipy import stats
35 +from sklearn.metrics import r2_score
36 +
37 +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
38 +from wp3 import config, design # noqa: E402
39 +
40 +GEO = config.RESULTS_DIR / "geo"
41 +PRED = GEO / "predictions"
42 +OUT = config.RESULTS_DIR / "geo_inference"
43 +DESIGNS = ["random", "size_matched", "ten_state", "gkf_state", "gkf_zip3", "gkf_zip5", "loso"]
44 +COMPARATORS = [("xgb", "full", "pred_nested"), ("xgb", "full", "pred_fixed1000"), ("xgb", "full", "pred_fixed100"),
45 + ("xgb", "plus_missing", "pred_nested"), ("lgb", "full", "pred_nested")]
46 +
47 +
48 +def load_pred(dsg, model, variant, seed):
49 + f = PRED / f"{dsg}_{model}_{variant}_seed{seed}.parquet"
50 + return pd.read_parquet(f) if f.exists() else None
51 +
52 +
53 +def seeds_available(dsg, model, variant):
54 + return sorted(int(p.stem.split("seed")[-1]) for p in PRED.glob(f"{dsg}_{model}_{variant}_seed*.parquet"))
55 +
56 +
57 +def within_r2(y, p, g):
58 + gs = pd.Series(g)
59 + yd = y - pd.Series(y).groupby(gs.values).transform("mean").values
60 + pdm = p - pd.Series(p).groupby(gs.values).transform("mean").values
61 + return r2_score(yd, pdm)
62 +
63 +
64 +def step_gap(info, draws):
65 + rng = np.random.RandomState(config.RANDOM_STATE)
66 + rows = []
67 + for dsg in DESIGNS:
68 + ols = load_pred(dsg, "ols", "full", 42)
69 + if ols is None:
70 + continue
71 + ols = ols.set_index("zpid")
72 + y = ols["y"].values
73 + z3 = info.loc[ols.index, "zip3"].values
74 + codes, uniq = pd.factorize(z3)
75 + members = [np.where(codes == g)[0] for g in range(len(uniq))]
76 + boot_idx = [np.concatenate([members[g] for g in rng.randint(0, len(uniq), len(uniq))]) for _ in range(draws)]
77 + base = {"design": dsg, "n_test": len(y), "n_zip3": len(uniq), "r2_ols": r2_score(y, ols["pred"].values)}
78 + for model, variant, col in COMPARATORS:
79 + seeds = seeds_available(dsg, model, variant)
80 + if not seeds:
81 + continue
82 + r2s, gaps = [], []
83 + for s in seeds:
84 + p = load_pred(dsg, model, variant, s).set_index("zpid").loc[ols.index]
85 + r2s.append(r2_score(y, p[col].values)); gaps.append(base["r2_ols"] - r2s[-1])
86 + # ZIP3-cluster bootstrap of the gap for the first seed (paired: same resampled listings)
87 + p = load_pred(dsg, model, variant, seeds[0]).set_index("zpid").loc[ols.index]
88 + po, pm = ols["pred"].values, p[col].values
89 + bg = np.array([r2_score(y[b], po[b]) - r2_score(y[b], pm[b]) for b in boot_idx])
90 + rows.append({**base, "comparator": f"{model}_{variant}_{col.replace('pred_', '')}", "n_seeds": len(seeds),
91 + "r2_comparator_mean": np.mean(r2s), "r2_comparator_sd": np.std(r2s, ddof=1) if len(r2s) > 1 else np.nan,
92 + "r2_comparator_min": np.min(r2s), "r2_comparator_max": np.max(r2s),
93 + "gap_ols_minus_comparator_mean": np.mean(gaps), "gap_sd_over_seeds": np.std(gaps, ddof=1) if len(gaps) > 1 else np.nan,
94 + "gap_min_over_seeds": np.min(gaps), "gap_max_over_seeds": np.max(gaps),
95 + "gap_seed42": gaps[0], "gap_boot_se_zip3": bg.std(ddof=1), "gap_ci_low_zip3": np.percentile(bg, 2.5),
96 + "gap_ci_high_zip3": np.percentile(bg, 97.5), "n_draws": draws})
97 + print(f" {dsg:12s} {rows[-1]['comparator']:28s} OLS {base['r2_ols']:.4f} comp {np.mean(r2s):.4f} "
98 + f"gap {np.mean(gaps):+.4f} [{rows[-1]['gap_ci_low_zip3']:+.4f}, {rows[-1]['gap_ci_high_zip3']:+.4f}]", flush=True)
99 + out = pd.DataFrame(rows); out.to_csv(OUT / "gap_bootstrap.csv", index=False)
100 + return out
101 +
102 +
103 +def step_loso(info):
104 + ols = load_pred("loso", "ols", "full", 42)
105 + if ols is None:
106 + return
107 + ols = ols.set_index("zpid"); y = ols["y"].values; st = info.loc[ols.index, "address_state"].values
108 + per = []
109 + for s_ in sorted(np.unique(st)):
110 + m = st == s_
111 + row = {"state": s_, "n_test": int(m.sum()), "r2_ols": r2_score(y[m], ols["pred"].values[m]),
112 + "r2_within_ols": r2_score(y[m] - y[m].mean(), ols["pred"].values[m] - ols["pred"].values[m].mean()),
113 + "mean_error_ols": float(np.mean(ols["pred"].values[m] - y[m]))}
114 + for model, variant, col in COMPARATORS:
115 + seeds = seeds_available("loso", model, variant)
116 + if not seeds:
117 + continue
118 + vals, within, me = [], [], []
119 + for s in seeds:
120 + p = load_pred("loso", model, variant, s).set_index("zpid").loc[ols.index][col].values
121 + vals.append(r2_score(y[m], p[m])); within.append(r2_score(y[m] - y[m].mean(), p[m] - p[m].mean())); me.append(np.mean(p[m] - y[m]))
122 + k = f"{model}_{variant}_{col.replace('pred_', '')}"
123 + row[f"r2_{k}"] = np.mean(vals); row[f"r2_{k}_sd"] = np.std(vals, ddof=1) if len(vals) > 1 else np.nan
124 + row[f"r2_within_{k}"] = np.mean(within); row[f"mean_error_{k}"] = np.mean(me)
125 + row[f"diff_ols_minus_{k}"] = row["r2_ols"] - np.mean(vals)
126 + per.append(row)
127 + per = pd.DataFrame(per); per.to_csv(OUT / "loso_per_state.csv", index=False)
128 + tests, summ = [], []
129 + cols = [c for c in per.columns if c.startswith("r2_") and not c.endswith("_sd") and "within" not in c]
130 + for c in cols:
131 + v = per[c].values; w = per["n_test"].values
132 + summ.append({"model": c.replace("r2_", ""), "median": np.median(v), "listing_weighted_mean": np.average(v, weights=w),
133 + "unweighted_mean": v.mean(), "p10": np.percentile(v, 10), "p90": np.percentile(v, 90), "min": v.min(), "max": v.max(),
134 + "n_states_gt_05": int((v > 0.5).sum()), "n_negative": int((v < 0).sum()),
135 + "median_abs_mean_error": float(np.median(np.abs(per[c.replace("r2_", "mean_error_")].values))) if c.replace("r2_", "mean_error_") in per else np.nan})
136 + if c != "r2_ols":
137 + d = per["r2_ols"].values - v
138 + nz = d[d != 0]
139 + tests.append({"comparison": f"OLS minus {c.replace('r2_', '')}", "n_states": len(d), "mean_diff": d.mean(), "median_diff": np.median(d),
140 + "n_states_ols_better": int((d > 0).sum()), "sign_test_p": float(stats.binomtest(int((nz > 0).sum()), len(nz)).pvalue),
141 + "wilcoxon_p": float(stats.wilcoxon(d).pvalue), "paired_t_p": float(stats.ttest_1samp(d, 0).pvalue),
142 + "listing_weighted_mean_diff": float(np.average(d, weights=w))})
143 + pd.DataFrame(summ).to_csv(OUT / "loso_summaries.csv", index=False)
144 + pd.DataFrame(tests).to_csv(OUT / "loso_paired_tests.csv", index=False)
145 + print(pd.DataFrame(summ).round(3).to_string(index=False)); print(pd.DataFrame(tests).round(4).to_string(index=False))
146 +
147 +
148 +def step_decomposition(info):
149 + rows = []
150 + for dsg in DESIGNS:
151 + for model, variant, col in [("ols", "full", "pred"), ("ols", "no_region", "pred"), ("ols_state_fe", "no_region", "pred")] + COMPARATORS:
152 + seeds = seeds_available(dsg, model, variant)
153 + if not seeds:
154 + continue
155 + p = load_pred(dsg, model, variant, seeds[0])
156 + y = p["y"].values; pr = p[col].values
157 + st = info.loc[p["zpid"].values, "address_state"].values
158 + groups = {"state": st}
159 + if dsg == "gkf_zip3":
160 + groups["zip3"] = info.loc[p["zpid"].values, "zip3"].values
161 + if dsg == "gkf_zip5":
162 + groups["zip5"] = info.loc[p["zpid"].values, "zip5"].values
163 + row = {"design": dsg, "model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "seed": seeds[0], "n": len(y), "r2": r2_score(y, pr)}
164 + e = pr - y; sst = ((y - y.mean()) ** 2).sum()
165 + for gname, g in groups.items():
166 + eb = pd.Series(e).groupby(pd.Series(g).values).transform("mean").values
167 + yd = y - pd.Series(y).groupby(pd.Series(g).values).transform("mean").values
168 + row[f"r2_within_{gname}_oracle"] = within_r2(y, pr, g)
169 + row[f"sse_within_share_{gname}"] = ((e - eb) ** 2).sum() / sst
170 + row[f"sse_between_share_{gname}"] = (eb ** 2).sum() / sst
171 + row[f"between_share_of_variance_{gname}"] = 1 - yd.var() / y.var()
172 + rows.append(row)
173 + out = pd.DataFrame(rows); out.to_csv(OUT / "sse_decomposition.csv", index=False)
174 + print(out[["design", "model", "r2", "r2_within_state_oracle", "sse_within_share_state", "sse_between_share_state"]].round(4).to_string(index=False))
175 +
176 +
177 +def step_feasible(info, ks=(5, 20, 50, 100, 200), reps=100):
178 + rng = np.random.RandomState(config.RANDOM_STATE)
179 + rows = []
180 + for dsg in ["loso", "ten_state"]:
181 + for model, variant, col in [("ols", "full", "pred"), ("xgb", "full", "pred_nested"), ("xgb", "full", "pred_fixed1000")]:
182 + seeds = seeds_available(dsg, model, variant)
183 + if not seeds:
184 + continue
185 + p = load_pred(dsg, model, variant, seeds[0]); y = p["y"].values; pr = p[col].values
186 + st = info.loc[p["zpid"].values, "address_state"].values
187 + states = sorted(np.unique(st))
188 + for k in ks:
189 + r2_pooled, r2_states = [], {s_: [] for s_ in states}
190 + for rep in range(reps):
191 + adj = pr.copy()
192 + for s_ in states:
193 + m = np.where(st == s_)[0]
194 + if len(m) <= k:
195 + continue
196 + cal = rng.choice(m, k, replace=False)
197 + level = np.mean(y[cal] - pr[cal])
198 + adj[m] = pr[m] + level
199 + rest = np.setdiff1d(m, cal)
200 + r2_states[s_].append(r2_score(y[rest], adj[rest]))
201 + r2_pooled.append(r2_score(y, adj))
202 + # oracle = all listings used for the level
203 + adj_o = pr + (pd.Series(y - pr).groupby(st).transform("mean").values)
204 + rows.append({"design": dsg, "model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "k_calibration": k, "reps": reps,
205 + "r2_raw": r2_score(y, pr), "r2_oracle_level": r2_score(y, adj_o),
206 + "r2_feasible_median": np.median(r2_pooled), "r2_feasible_p25": np.percentile(r2_pooled, 25), "r2_feasible_p75": np.percentile(r2_pooled, 75),
207 + "r2_feasible_mean": np.mean(r2_pooled),
208 + "median_state_r2_feasible": np.median([np.median(v) for v in r2_states.values() if v])})
209 + print(f" {dsg} {rows[-1]['model']} k={k}: raw {rows[-1]['r2_raw']:.4f} oracle {rows[-1]['r2_oracle_level']:.4f} feasible {rows[-1]['r2_feasible_median']:.4f}", flush=True)
210 + pd.DataFrame(rows).to_csv(OUT / "feasible_level_correction.csv", index=False)
211 +
212 +
213 +def step_two_by_two(info):
214 + met = pd.read_csv(GEO / "metrics.csv")
215 + pooled = pd.read_csv(GEO / "metrics_pooled.csv") if (GEO / "metrics_pooled.csv").exists() else None
216 + rows = []
217 + for dsg in ["ten_state", "loso", "gkf_state"]:
218 + for variant in ["full", "no_region", "no_school", "no_school_no_region", "plus_missing"]:
219 + for rule in ["nested", "fixed100", "fixed1000"]:
220 + if dsg == "ten_state":
221 + sub = met[(met.design == dsg) & (met.model == "xgb") & (met.variant == variant) & (met.rounds_rule == rule)]
222 + if sub.empty:
223 + continue
224 + rows.append({"design": dsg, "variant": variant, "rounds_rule": rule, "region_dummies": "no_region" not in variant,
225 + "school_vars": "no_school" not in variant, "r2_mean": sub.r2.mean(), "r2_sd": sub.r2.std(ddof=1) if len(sub) > 1 else np.nan,
226 + "n_seeds": len(sub), "r2_within_state_oracle": sub.r2_within_state.mean(), "mdape": sub.mdape.mean(),
227 + "selected_rounds": sub.selected_rounds.iloc[0] if "selected_rounds" in sub else np.nan})
228 + else:
229 + sub = pooled[(pooled.design == dsg) & (pooled.model == "xgb") & (pooled.variant == variant) & (pooled.rounds_rule == rule)] if pooled is not None else pd.DataFrame()
230 + if sub.empty:
231 + continue
232 + per = met[(met.design == dsg) & (met.model == "xgb") & (met.variant == variant) & (met.rounds_rule == rule)]
233 + rows.append({"design": dsg, "variant": variant, "rounds_rule": rule, "region_dummies": "no_region" not in variant,
234 + "school_vars": "no_school" not in variant, "r2_mean": sub.r2.mean(), "r2_sd": sub.r2.std(ddof=1) if len(sub) > 1 else np.nan,
235 + "n_seeds": len(sub), "r2_within_state_oracle": sub.r2_within_state.mean(), "mdape": sub.mdape.mean(),
236 + "median_fold_r2": per.groupby("seed").r2.median().mean()})
237 + out = pd.DataFrame(rows); out.to_csv(OUT / "two_by_two_region_rounds.csv", index=False)
238 + print(out.round(4).to_string(index=False))
239 +
240 +
241 +def step_duplicates(info):
242 + dup = info["dup_group_id"].values
243 + rows = []
244 + for model, variant, col in [("ols", "full", "pred")] + COMPARATORS[:2]:
245 + p = load_pred("random", model, variant, 42)
246 + if p is None:
247 + continue
248 + y = p["y"].values; pr = p[col].values
249 + d = info.loc[p["zpid"].values, "dup_group_id"].values
250 + # a test row has a duplicate in train iff its group has a member not in the test set
251 + test_z = set(p["zpid"].values)
252 + grp_all = info[info["dup_group_id"] >= 0]
253 + in_train = grp_all[~grp_all.index.isin(test_z)].groupby("dup_group_id").size()
254 + has_train_dup = (d >= 0) & pd.Series(d).map(in_train).fillna(0).values.astype(bool)
255 + rows.append({"model": f"{model}_{variant}_{col.replace('pred_', '').replace('pred', 'ols')}", "n_test": len(y),
256 + "n_test_with_train_duplicate": int(has_train_dup.sum()), "r2_all": r2_score(y, pr),
257 + "r2_without_duplicates": r2_score(y[~has_train_dup], pr[~has_train_dup]),
258 + "r2_duplicates_only": r2_score(y[has_train_dup], pr[has_train_dup]) if has_train_dup.sum() > 10 else np.nan})
259 + pd.DataFrame(rows).to_csv(OUT / "duplicates_random_split_v2.csv", index=False); print(pd.DataFrame(rows).round(4).to_string(index=False))
260 +
261 +
262 +def step_range(gap):
263 + if gap is None or gap.empty:
264 + return
265 + sub = gap[gap.comparator.isin(["xgb_full_nested", "xgb_full_fixed1000", "xgb_full_fixed100", "lgb_full_nested"])]
266 + rows = []
267 + for comp, g in sub.groupby("comparator"):
268 + rows.append({"comparator": comp, "designs": " ".join(g.design), "gap_min_design": g.gap_ols_minus_comparator_mean.min(),
269 + "gap_max_design": g.gap_ols_minus_comparator_mean.max(),
270 + "design_of_min": g.loc[g.gap_ols_minus_comparator_mean.idxmin(), "design"],
271 + "design_of_max": g.loc[g.gap_ols_minus_comparator_mean.idxmax(), "design"],
272 + "min_ci_low": g.gap_ci_low_zip3.min(), "max_ci_high": g.gap_ci_high_zip3.max(),
273 + "n_designs_ci_excludes_zero": int(((g.gap_ci_low_zip3 > 0) | (g.gap_ci_high_zip3 < 0)).sum()), "n_designs": len(g)})
274 + pd.DataFrame(rows).to_csv(OUT / "gap_range.csv", index=False); print(pd.DataFrame(rows).round(4).to_string(index=False))
275 +
276 +
277 +def main():
278 + ap = argparse.ArgumentParser(description=__doc__)
279 + ap.add_argument("--draws", type=int, default=1000)
280 + ap.add_argument("--steps", nargs="+", default=["gap", "loso", "decomposition", "feasible", "two_by_two", "duplicates", "range"])
281 + args = ap.parse_args()
282 + OUT.mkdir(parents=True, exist_ok=True)
283 + t0 = time.time()
284 + info = pd.read_parquet(design.DESIGN_V2, columns=["zpid", "address_state", "zip3", "zip5", "dup_group_id"]).set_index("zpid")
285 + gap = None
286 + for s in args.steps:
287 + print(f"== {s} ({time.time() - t0:.0f}s)", flush=True)
288 + if s == "gap":
289 + gap = step_gap(info, args.draws)
290 + elif s == "range":
291 + step_range(gap if gap is not None else (pd.read_csv(OUT / "gap_bootstrap.csv") if (OUT / "gap_bootstrap.csv").exists() else None))
292 + else:
293 + {"loso": step_loso, "decomposition": step_decomposition, "feasible": step_feasible, "two_by_two": step_two_by_two,
294 + "duplicates": step_duplicates}[s](info)
295 + print(f"Done in {time.time() - t0:.0f}s")
296 +
297 +
298 +if __name__ == "__main__":
299 + main()
300
Menu

Theme
↑↓ navigate ↵ open esc close / or ⌘K anywhere