Python 空间自相关实战:70 城 20 年房价数据告诉你涨跌有没有规律
- 2026-06-28 10:51:27
Python 空间自相关实战:70 城 20 年房价数据告诉你涨跌有没有规律
累计涨幅 vs 波动率散点图 I = n ∑ i ∑ j w i j ⋅ ∑ i ∑ j w i j ( x i − x ˉ ) ( x j − x ˉ ) ∑ i ( x i − x ˉ ) 2
涨跌幅 vs 波动率的 Moran 散点图对比
这是一个系列的第一篇。本篇覆盖数据准备、探索性分析和全局空间自相关（Moran’s I），下篇讲局部分析（LISA）和地图可视化。
1. 你要回答的问题
拿到 70 个城市 20 年的房价数据，第一个直觉是："涨得多的城市是不是挨在一起？"这就是空间自相关——衡量"邻居之间是否相似"的统计方法。
本文用 Python 复现完整流程，最终发现一个反直觉的结果：涨跌幅在空间上没有聚集，但波动率有。
2. 环境准备
安装 uv 之后创建一个项目。
uv的安装办法[1]
uv init housing-spatial && cd housing-spatialuv add pandas geopandas numpy matplotlib libpysal esda adjustText shapelypandas | |
geopandas | |
libpysal | |
esda | |
matplotlib | |
shapely |
3. Step 1：数据获取与理解
数据来源：hugohe3/70cityprice[2]，国家统计局 70 城商品住宅销售价格指数，月度数据，2006-2026 年。
import pandas as pddf = pd.read_csv('70cityprice.csv')print(f'数据量: {len(df)} 行, {df["CITY"].nunique()} 个城市')print(f'时间范围: {df["DATE"].min()} ~ {df["DATE"].max()}')数据量: 47040 行, 70 个城市时间范围: 2006-01-01 ~ 2026-04-01数据长这样（以北京 2006 年前几个月为例）：
| 环比 | 101.1 | ||||
| 环比 | 101.3 |
每个城市每个月有 2-3 行，按 FixedBase 区分：
| 环比 | |||
为什么只用环比？ 环比是月度价格变化的直接度量。同比消除了季节性但也消除了月度波动信息。定基比数据不完整。我们要算的"波动率"就是月度涨跌的幅度，用环比才对。
CommodityHouseIDX 是新建商品住宅指数，SecondHandIDX 是二手住宅指数。本教程只用 CommodityHouseIDX。
4. Step 2：从环比指数到累计涨幅和波动率
df['DATE'] = pd.to_datetime(df['DATE'])df_mom = df[df['FixedBase'] == '环比'].copy()df_mom['index_value'] = pd.to_numeric(df_mom['CommodityHouseIDX'], errors='coerce')stats = []for city in df_mom['CITY'].unique(): c = df_mom[df_mom['CITY'] == city].sort_values('DATE') cum_return = (c['index_value'] / 100).prod() * 100 - 100# 累计涨幅 volatility = c['index_value'].std() # 波动率（标准差） stats.append({'CITY': city, 'cum_return': cum_return, 'volatility': volatility})stats_df = pd.DataFrame(stats)print(stats_df.sort_values('cum_return', ascending=False).head(5)) CITY cum_return volatility0 北京 206.2% 1.07%1 上海 197.4% 1.02%2 西安 190.3% 0.82%3 三亚 186.8% 2.04%4 海口 184.6% 1.82%累计涨幅公式：
直觉：如果每个月都涨 1%（环比=101），244 个月后累计涨
4.1 探索性分析：波动率 vs 涨幅
import matplotlib.pyplot as pltfig, ax = plt.subplots(figsize=(10, 7))ax.scatter(stats_df['volatility'], stats_df['cum_return'], s=50, alpha=0.6)# 标注典型城市for _, r in stats_df.iterrows():if r['CITY'] in ['北京', '上海', '深圳', '三亚', '海口', '温州', '成都', '北海']: ax.annotate(r['CITY'], (r['volatility'], r['cum_return']), fontsize=9, fontweight='bold', xytext=(5, 5), textcoords='offset points')r = stats_df['volatility'].corr(stats_df['cum_return'])ax.text(0.05, 0.95, f'r = {r:.3f}', transform=ax.transAxes, fontsize=11, va='top')ax.set_xlabel('波动率（月度环比指数标准差）')ax.set_ylabel('20年累计涨幅（%）')ax.set_title('70城房价：波动率 vs 累计涨幅')ax.grid(True, alpha=0.2)plt.savefig('vol_vs_return.png', dpi=150, bbox_inches='tight')波动率和累计涨幅相关系数 r = 0.55——高波动城市倾向于涨得多，但不完全重合。温州是高波动低涨幅（民间借贷危机），成都是低波动高涨幅（稳步上涨）。
5. Step 3：把城市数据变成空间数据
Moran’s I 需要知道每个城市的经纬度，才能算谁跟谁是邻居。CSV 只有城市名，没有坐标。
70 城的坐标需要手动收集或用 API 地理编码。以下是完整字典：
CITY_COORDS = {"三亚": (109.51, 18.25), "上海": (121.48, 31.22), "丹东": (124.37, 40.13),"乌鲁木齐": (87.68, 43.77), "九江": (115.99, 29.71), "兰州": (103.73, 36.03),"包头": (109.84, 40.65), "北京": (116.46, 39.92), "北海": (109.12, 21.48),"南京": (118.78, 32.04), "南充": (106.08, 30.83), "南昌": (115.89, 28.68),"南宁": (108.33, 22.84), "厦门": (118.10, 24.46), "合肥": (117.27, 31.86),"哈尔滨": (126.63, 45.75), "唐山": (118.02, 39.63), "吉林": (126.57, 43.87),"呼和浩特": (111.65, 40.82), "咸阳": (108.71, 34.33), "大连": (121.62, 38.92),"天津": (117.10, 39.10), "太原": (112.55, 37.87), "宁波": (121.56, 29.87),"安庆": (117.03, 30.52), "宜昌": (111.30, 30.70), "宝鸡": (107.14, 34.37),"平顶山": (113.19, 33.77), "广州": (113.23, 23.16), "徐州": (117.20, 34.26),"成都": (104.06, 30.67), "扬州": (119.42, 32.39), "承德": (117.93, 40.97),"拉萨": (91.11, 29.97), "无锡": (120.29, 31.59), "昆明": (102.73, 25.04),"杭州": (120.19, 30.26), "桂林": (110.29, 25.27), "武汉": (114.31, 30.52),"泸州": (105.39, 28.91), "洛阳": (112.44, 34.63), "济南": (117.00, 36.65),"济宁": (116.59, 35.41), "海口": (110.35, 20.02), "深圳": (114.07, 22.62),"温州": (120.65, 28.01), "湖州": (120.09, 30.87), "滁州": (118.32, 32.25),"烟台": (121.39, 37.52), "牡丹江": (129.58, 44.55), "石家庄": (114.48, 38.03),"福州": (119.30, 26.08), "秦皇岛": (119.57, 39.95), "秦皇岛": (119.57, 39.95),"绍兴": (120.58, 30.00), "绵阳": (104.73, 31.47), "芜湖": (118.38, 31.33),"蚌埠": (117.34, 32.92), "衡水": (115.72, 37.72), "赣州": (114.94, 25.85),"遵义": (106.90, 27.73), "邢台": (114.50, 37.07), "郑州": (113.65, 34.76),"重庆": (106.54, 29.59), "金华": (119.64, 29.12), "银川": (106.27, 38.47),"锦州": (121.15, 41.13), "长春": (125.35, 43.88), "长沙": (112.94, 28.23),"青岛": (120.33, 36.07), "韶关": (113.62, 24.84), "湛江": (110.36, 21.27),"大理": (100.19, 25.69), "泉州": (118.58, 24.93), "岳阳": (113.09, 29.37),"常德": (111.69, 29.05), "惠州": (114.42, 23.09), "沈阳": (123.38, 41.80),"泸州": (105.39, 28.91), "贵阳": (106.71, 26.57), "赣州": (114.94, 25.85),"西安": (108.95, 34.27), "襄阳": (112.14, 32.04), "西宁": (101.74, 36.56),}stats_df['lng'] = stats_df['CITY'].map(lambda c: CITY_COORDS.get(c, (None, None))[0])stats_df['lat'] = stats_df['CITY'].map(lambda c: CITY_COORDS.get(c, (None, None))[1])missing = stats_df[stats_df['lng'].isna()]['CITY'].tolist()if missing:print(f'⚠ 缺少坐标的城市: {missing}') stats_df = stats_df.dropna(subset=['lng', 'lat'])转成 GeoDataFrame：
import geopandas as gpdfrom shapely.geometry import Pointgeometry = [Point(xy) for xy inzip(stats_df['lng'], stats_df['lat'])]gdf = gpd.GeoDataFrame(stats_df, geometry=geometry, crs='EPSG:4326')gdf.to_file('housing_data.geojson', driver='GeoJSON')print(f'GeoDataFrame: {len(gdf)} 个城市, CRS = {gdf.crs}')CRS = EPSG:4326 是 WGS84 坐标系（GPS 用的），经度 [-180, 180]，纬度 [-90, 90]。空间分析必须明确坐标系，否则距离计算会出错。
6. Step 4：空间自相关是什么？为什么要算它？
核心问题：北京涨了 206%，天津涨了多少？如果天津也涨了很多，说明"涨"在空间上传导了；如果天津没怎么涨，说明各城市独立。
空间自相关就是量化这种"邻居之间是否相似"的统计量。
6.1 Moran’s I 的公式
| 协变 |
直觉：涨得多的城市（
- I > 0
：高的挨着高的（正聚集） - I ≈ 0
：随机 - I < 0
：高的挨着低的（棋盘格局） 理论范围 [-1, 1]，实际很少超过 0.5
7. Step 5：构建空间权重矩阵
70 城是点数据，分布不均匀（东部密、西部稀），用 KNN（k 近邻）：
import libpysal as lpsw = lps.weights.KNN.from_dataframe(gdf, k=5)w.transform = 'r'# 行标准化# 看看北京的邻居是谁print(f'北京的 5 个邻居: {gdf.loc[list(w.neighbors[gdf[gdf["CITY"]=="北京"].index[0]]), "CITY"].values}')为什么 KNN 而不是距离阈值？ 70 城分布极不均匀。用 500km 阈值，北京可能有 10 个邻居，乌鲁木齐可能 0 个。KNN 保证每个城市恰好有 k 个邻居。
为什么 k=5？ 经验值：样本量的 5-10%。70 × 7% ≈ 5。后面会做 k=3/5/7 的敏感性检验。
行标准化：不标准化的话，每个邻居权重=1，东部城市邻居密集、西部稀疏，尺度不一致。标准化后每个城市的邻居权重之和 = 1。
8. Step 6：计算 Moran’s I
from esda.moran import Morany_cum = gdf['cum_return'].valuesy_vol = gdf['volatility'].valuesmoran_cum = Moran(y_cum, w)moran_vol = Moran(y_vol, w)print(f"累计涨幅: I = {moran_cum.I:.4f}, p = {moran_cum.p_sim:.4f}")print(f"波动率: I = {moran_vol.I:.4f}, p = {moran_vol.p_sim:.4f}")累计涨幅: I = -0.0065, p = 0.4420波动率: I = 0.3016, p = 0.00108.1 p 值是怎么算的？排列检验
p_sim = 0.001 是蒙特卡洛模拟：
把 70 个城市的值随机打乱（空间位置不变） 算一次 Moran’s I 重复 999 次 看真实 I 在 999 次中排第几
# 查看排列检验细节print(f'排列次数: {moran_vol.permutations}')print(f'模拟分布均值: {moran_vol.EI:.4f}')print(f'z 值: {moran_vol.z_sim:.4f}') # >1.96 即 p<0.059. Step 7：画 Moran 散点图对比
import numpy as npfig, axes = plt.subplots(1, 2, figsize=(16, 7))for ax, y, title in [(axes[0], y_cum, '累计涨幅'), (axes[1], y_vol, '波动率')]: y_std = (y - y.mean()) / y.std() y_lag = lps.weights.lag_spatial(w, y_std)# 判断象限 types = []for zs, zl inzip(y_std, y_lag):if zs >= 0and zl >= 0: types.append('HH')elif zs < 0and zl < 0: types.append('LL')elif zs < 0and zl >= 0: types.append('LH')else: types.append('HL') colors = {'HH': '#e15759', 'LL': '#4e79a7', 'LH': '#76b7b2', 'HL': '#f28e2b'}for t in ['HH', 'LL', 'LH', 'HL']: mask = [i for i, tt inenumerate(types) if tt == t] ax.scatter(y_std[mask], y_lag[mask], c=colors[t], label=t, s=50, alpha=0.7)# 回归线 z_coef = np.polyfit(y_std, y_lag, 1) x_line = np.linspace(y_std.min(), y_std.max(), 100) ax.plot(x_line, np.poly1d(z_coef)(x_line), 'r--', lw=2) ax.axhline(0, color='black', lw=0.5) ax.axvline(0, color='black', lw=0.5) moran = Moran(y, w) sig = '显著'if moran.p_sim < 0.05else'不显著' ax.set_title(f'{title} (I={moran.I:.3f}, p={moran.p_sim:.3f}, {sig})', fontweight='bold') ax.set_xlabel('标准化值') ax.set_ylabel('空间滞后') ax.legend() ax.grid(True, alpha=0.2)# 标注典型城市for city_name in ['三亚', '海口', '北京', '温州', '成都']: idx = gdf[gdf['CITY'] == city_name].indexiflen(idx) > 0: i = idx[0] ax.annotate(city_name, (y_std[i], y_lag[i]), fontsize=8, fontweight='bold', xytext=(6, 6), textcoords='offset points')plt.suptitle('Moran散点图对比：涨跌幅 vs 波动率', fontsize=15, fontweight='bold')plt.tight_layout()plt.savefig('moran_comparison.png', dpi=150, bbox_inches='tight')左图（累计涨幅）：点是一团乱麻，回归线接近水平。右图（波动率）：点明显向右上倾斜，HH 象限有三亚、海口、深圳扎堆。
10. Step 8：敏感性检验
for k in [3, 5, 7]: w_k = lps.weights.KNN.from_dataframe(gdf, k=k) w_k.transform = 'r' m_cum = Moran(y_cum, w_k) m_vol = Moran(y_vol, w_k)print(f'k={k}: 累计涨幅 I={m_cum.I:.4f} p={m_cum.p_sim:.4f} | 波动率 I={m_vol.I:.4f} p={m_vol.p_sim:.4f}')k=3: 累计涨幅 I=-0.0039 p=0.4830 | 波动率 I= 0.4119 p=0.0010k=5: 累计涨幅 I=-0.0065 p=0.4420 | 波动率 I= 0.3016 p=0.0020k=7: 累计涨幅 I=-0.0076 p=0.4700 | 波动率 I= 0.3266 p=0.0010结论稳健：无论 k=3/5/7，波动率的 Moran’s I 都在 0.30-0.41 之间，p 值始终 < 0.002。
11. 本篇小结
一句话：买房不用看邻居涨了多少，但如果你的城市波动率高（比如海南北部湾），周边城市大概率也高。
12. 下篇预告
LISA 局部空间自相关——Moran’s I 只告诉你"有没有聚集"，不告诉你"谁跟谁聚集"。LISA 能定位具体哪些城市在扎堆。三亚、海口、北海、湛江为什么形成全国唯一的波动率热点？
参考链接
[1] https://docs.astral.sh/uv/getting-started/installation/
[2] https://github.com/hugohe3/70cityprice
[3] https://github.com/renhai-lab/Urban-Spatial-Data-Analysis-Notebook/tree/master/projects/housing-analysis
本文来自网友投稿或网络内容，如有侵犯您的权益请联系我们删除，联系邮箱：wyl860211@qq.com 。
