Python 空间自相关实战：70 城 20 年房价数据告诉你涨跌有没有规律

renhai

发布于 2026-06-23 17:42:34

1260

这是一个系列的第一篇。本篇覆盖数据准备、探索性分析和全局空间自相关（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 shapely

库	干什么用
pandas	数据处理
geopandas	空间数据（带坐标信息的 DataFrame）
libpysal	空间权重矩阵（定义"谁是谁的邻居"）
esda	空间自相关统计量（Moran’s I、LISA）
matplotlib	画图
shapely	几何对象（Point）

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 年前几个月为例）：

DATE	ADCODE	CITY	FixedBase	CommodityHouseIDX	SecondHandIDX
2006-01	110100	北京	同比	108.4	110.2
2006-01	110100	北京	环比	101.1	101.4
2006-02	110100	北京	同比	107.3	110.1
2006-02	110100	北京	环比	101.3	102.3

每个城市每个月有 2-3 行，按 FixedBase 区分：

FixedBase	含义	基准	本教程是否用
环比	本月 vs 上月	上月 = 100	✅ 用这个
同比	本月 vs 去年同月	去年同月 = 100	❌
定基比	本月 vs 固定基期	基期 = 100	❌（仅 2011 后有）

为什么只用环比？ 环比是月度价格变化的直接度量。同比消除了季节性但也消除了月度波动信息。定基比数据不完整。我们要算的"波动率"就是月度涨跌的幅度，用环比才对。

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%

累计涨幅公式：累计涨幅=(∏t=1244It100−1)×100%累计涨幅=(∏t=1244100It−1)×100%

直觉：如果每个月都涨 1%（环比=101），244 个月后累计涨 1.01244−1=1066%1.01244−1=1066%。实际北京的月度环比均值约 100.47，连乘后累计涨 206%。

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')

累计涨幅 vs 波动率散点图

波动率和累计涨幅相关系数 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=n∑i∑jwij⋅∑i∑jwij(xi−xˉ)(xj−xˉ)∑i(xi−xˉ)2I=∑i∑jwijn⋅∑i(xi−xˉ)2∑i∑jwij(xi−xˉ)(xj−xˉ)

符号	含义
xixi	城市 i 的值（比如累计涨幅）
xˉxˉ	所有城市的均值
wijwij	空间权重——城市 i 和 j 是否为邻居
(xi−xˉ)(xj−xˉ)(xi−xˉ)(xj−xˉ)	协变——i 和 j 的偏差是否同号

直觉：涨得多的城市（xi>xˉxi>xˉ）的邻居也涨得多（xj>xˉxj>xˉ），乘积为正 → Moran’s I > 0 → 正聚集。

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.0010

8.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.05

9. 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].index        iflen(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')

涨跌幅 vs 波动率的 Moran 散点图对比

左图（累计涨幅）：点是一团乱麻，回归线接近水平。右图（波动率）：点明显向右上倾斜，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。