第 7 章 · 02 数据清洗与合成控制法

第 7 章 · 02 数据清洗与合成控制法

本节摘要:本节精读 Smart Farmers 的数据层——cleanse data.py(约 363 行)和 country selection.py(约 308 行)。数据从 FAO(联合国粮农组织)数据库下载,解压后约 14GB,涵盖生产、产量、价格、人口、人均 GDP 等维度。我们看作者怎么把马来西亚从这堆巨量数据里捞出来,先按"用地占比 <1% 的剔除"把 61 类作物砍到 30 类,再用 mapping.csv 作物寿命表给多年生作物套上寿命约束。本节真正的亮点是合成控制法(Synthetic Control Method)——对 FAO 缺失价格的作物(姜、生菜、玉米、橙、甘蔗、烟草),作者用其他国家作控制组加权拟合出"合成马来西亚价格",灵感直接来自 Abadie & Gardeazabal (2003) 那篇用合成控制法评估巴斯克恐怖主义经济成本的经典论文。本节末尾还会点出工程层的两个坑:脚本对 global 变量的重度依赖、旧版 pandas append API 在新版会告警。

涉及源码:原项目 Smart Farmers project/cleanse data.py(约 363 行),Smart Farmers project/country selection.py(约 308 行),辅助文件 mapping.csv(作物寿命表)。

⚠️ 注意:本节代码是典型的 Jupyter Notebook 导出风格(.py 里残留大量 # In[1]: cell 标记),并且重度依赖 global 变量(beginyear / endyear / eco_coeff / D 等)。这意味着函数签名虽然干净,但脱离了脚本顶部的 global 声明就跑不动。读代码时务必自上而下,不要孤立看函数。

学习目标

阅读完本节,你应当能够:

1. 说清 FAO 数据库的五个核心维度(Production / Prices / LandUse / Population / Macro)和它们的 join key。
2. 解释"61 类作物 → 30 类"的两步筛选(用地 <1% 剔除 + 价格缺失剔除)。
3. 复述多年生作物寿命约束的公式 (ωL−1)/(ωL) 及 mapping.csv 的取数规则(Google 首页结果取均值)。
4. 讲清合成控制法在 Smart Farmers 里的角色转换——从"评估政策效果"变成"补缺失价格"。
5. 看懂 country selection.py 怎么用"出口占比"挑控制组国家。
6. 知道脚本里 global 满天飞和 malay_gdp.append 旧 API 两个工程坑。

金融/数学原理

FAO:一个 14GB 的农业数据黑洞

README 的 Data Specification 段开宗明义:

we can find most of the data with consistent format from one source, Food and Agriculture Organization. FAO is an organization under United Nations and it provides a database with absurd level of details. The database is about 14 gigabytes after decompression. The csv of international trade alone takes up 4 gigabytes.

14GB 是什么概念?本章之前所有的策略数据加起来(K 线、宏观数据、期货价格)不超过几百 MB。FAO 的"absurd level of details"包括:

• Production_Crops_E_All_Data:各作物各年的种植面积(Area harvested)与产量(Production,吨);
• Prices_E_All_Data:各作物各年的生产者价格(Producer Price,USD/吨);
• Inputs_LandUse_E_All_Data:各国土地利用(Cropland 耕地总面积);
• Population_E_All_Data:各国人口(后续作需求回归的自变量);
• Macro-Statistics_Key_Indicators_E_All_Data:各国人均 GDP(需求回归的另一个自变量)。

cleansing 的任务就是把这五张表按"年份 + 国家 + 作物"三个 key join 起来,拼出一张大宽表。代码顶部的 pd.read_csv 一口气读了五个 csv:

#cleanse data.py:104-117(节选)
prod       = pd.read_csv('Production_Crops_E_All_Data_(Normalized).csv', encoding='latin-1')
prix       = pd.read_csv('Prices_E_All_Data_(Normalized).csv',          encoding='latin-1')
land       = pd.read_csv('Inputs_LandUse_E_All_Data_(Normalized).csv',  encoding='latin-1')
population = pd.read_csv('Population_E_All_Data_(Normalized).csv',      encoding='latin-1')
gdp        = pd.read_csv('Macro-Statistics_Key_Indicators_E_All_Data_(Normalized).csv', encoding='latin-1')
​

注意 encoding='latin-1'——FAO 数据里有大量带重音的国家名(如 Côte d'Ivoire),用 utf-8 读会炸,latin-1 是兜底编码。

选马来西亚:不是因为它爱吃椰浆饭

README 解释为何挑马来西亚做样本国:

The target country is handpicked as Malaysia because I am a big fan of Nasi Lemak. Lol no, I am not. The actual justification is its low variety of crops. For Malaysia, we only need to examine 61 categories of crops but for Italy the number becomes 111.

真正的理由是作物种类少 → 模型维度低 → QP 容易解。马来西亚 61 类,意大利 111 类,差别近一倍。30 维 QP 对 CVXOPT 是小菜,111 维虽然也能解,但调参会更痛苦。这是一个很务实的"先简化再推广"的研究策略。

时间窗方面,脚本里 beginyear=2012, endyear=2019(注意 2012 是上一年的缓冲,真正回测从 2013 起,因为定价机制要用 t-1 的价格做基准)。2013-2018 共 6 年做回测,2019-2025 做前向预测(下一节详谈)。

作物筛选:61 → 30 的两刀

cleansing 的第一刀是用地剔除。malay_prod 里很多" subtotal "类(如"Cereals, Total""Fruit Primary"),它们只是细分类的汇总,会和细分类重复计算。check consistency.py 里有完整的 exclude 清单(29 项 subtotal):

#check consistency.py:76-104(节选)
exclude=['Areca nuts','Bastfibres, other','Cashew nuts, with shell',
         'Cereals, Total','Chillies and peppers, dry','Citrus Fruit, Total',
         'Cloves','Coarse Grain, Total','Coffee, green','Coir', ...]
​

cleanse data.py 里还有第二刀——按 Item Code 剔除用地 <1% 的边角作物:

#cleanse data.py:151-156
exclude=[813,236,809,1717, 512,782,656,149,667,809,
         813,689,698,702,723,217, 226,236,242,463,603,619,
         1720,1729,1731,1732,1735, 1738,1753,1804,1841,1814]
​

两刀砍完,61 类 → 30 类。这 30 类就是后续 QP 的 30 维决策变量。

多年生作物寿命表:mapping.csv

cleansing 的一个关键 join 是 malay_prod.merge(mapping, on=['Item Code','Item'], how='left'),把作物寿命 L 接上去。README 对 mapping.csv 的取数规则特别坦诚:

If you search the lifespan of an apple tree, you will end up with multiple answers from 60 to 150 years. Let's define a tidy structure for this challenge. We only collect the information of the first page of google organic results and take the mean of different numbers. After that, we create mapping.csv and upload to GitHub. You are welcome to raise an issue when you disagree with any number inside this file.

"Google 首页结果取均值"听起来很糙,但作者把它显式声明出来,反而比假装精确更诚实——读者可以提 issue 改数。这是学术写作的好习惯:承认数据的局限比掩饰它更有价值。

寿命 L 拿到后,cleansing 里这行算的就是第 01 节讲过的多年生下界系数:

#cleanse data.py:73
perennial = D[currentyear]['lifespan'].dropna().apply(
    lambda x: (x*eco_coeff - 1) / (x*eco_coeff)
)
​

其中 eco_coeff=0.7(农学家建议的 ω)。对于一年生作物(lifespan 为 NaN),用默认惯性系数 default_discount=0.8 填充。

合成控制法:从"评估政策"到"补缺失价格"

本节最精彩的部分。FAO 数据虽全,但有些作物的生产者价格在马来西亚是缺失的——姜(ginger)、生菜(lettuce)、玉米(maize)、橙(orange)、甘蔗(sugarcane)、烟草(tobacco)。怎么补?

传统办法有四种:前值填充、后值填充、算术平均、卡尔曼滤波。作者偏偏选了第五种——合成控制法(Synthetic Control Method)。这个方法原本是政治学/政策评估的利器,Abadie & Gardeazabal (2003) 用它评估巴斯克恐怖活动对当地 GDP 的冲击:找一组"未受恐怖活动影响"的控制组地区,加权拟合出一个"合成巴斯克",对比真实巴斯克与合成巴斯克的 GDP 差,就是恐怖活动的经济成本。

Smart Farmers 把这个方法角色翻转——原本的"控制组"在这里是"有价格数据的其他国家","处理组"是"缺失价格的马来西亚",所谓的"政策"就是"价格数据缺失"。逻辑链很美:

• 农产品市场不是纯国内的,跨国套利会让同种作物在各国的价格趋同(README 假设 No arbitrage,这里正是这个假设的兑现);
• 所以马来西亚的价格可以用其他国家价格的加权组合来近似;
• 权重由 OLS 拟合得到(在马来西亚有数据的年份上拟合,再预测缺失年份)。

cleanse data.py 末尾被注释掉的代码就是原始的合成控制拟合过程:

#cleanse data.py:333-362(原始 SCM 拟合,注释保留)
#create synthetic control unit
for ii in ss['Item'].unique():
    temp = prix[prix['Item']==ii][...].pivot(index='Year', columns='Area', values='Value')
    # 删掉任何一年为空的国家列(保留 Malaysia)
    for i in temp:
        if temp[i].isnull().any() and i != 'Malaysia':
            del temp[i]
    # 训练集:Malaysia 有价格的年份
    x = temp.loc[temp['Malaysia'].dropna().index]; del x['Malaysia']
    y = temp['Malaysia'].dropna()
    m = sm.OLS(y, x).fit()
    # 预测:Malaysia 缺价格的年份
    test = temp[temp['Malaysia'].isnull()]
    print(m.predict(test).tolist())
​

拟合结果作者直接硬编码进了脚本(L279-294):

#cleanse data.py:279-294
malay_ginger    = [1711.63, 1918.06, 1926.35, 1904.05]
malay_lettuce   = [808.91, 809.29, 880.07, 801.19]
malay_maize     = [176.03, 181.61, 185.21]
malay_orange    = [385.21, 303.63, 277.52, 376.93]
malay_sugarcane = [246.78]
malay_tobacco   = [4345.96, 4091.37, 3638.90, 3829.53, 3758.18, 3889.60]
​

这些数字就是用其他国家作控制组、OLS 加权拟合出的"合成马来西亚价格"。注意它们只覆盖有缺失的年份,有的作物缺 4 年(姜)、有的只缺 1 年(甘蔗)、烟草缺 6 年(2013-2018 全缺)。这一步做完,所有 30 类作物的价格时间序列就齐全了。

country selection.py:用出口占比筛控制组

合成控制法的关键是"挑哪些国家当控制组"。country selection.py 就是为这件事服务的——它计算各国"出口量占国内产量"的比例(出口依赖度),按这个指标排序。

#country selection.py:131-166(节选)
export = trade.groupby(['Area','Year']).sum()       # 各国年出口总量
supply = prod.groupby(['Area','Year']).sum()        # 各国年产量
pourcent['value'] = np.divide(export['Value'], supply['Value'])   # 出口/产量
mean_pourcent = pourcent[['Area','value']].groupby('Area').mean() # 历史均值
mean_pourcent = mean_pourcent.sort_values('value')  # 排序
​

为什么用出口占比?因为合成控制法需要"和目标国经济结构相近"的控制组。出口占比低(自给自足型)的国家,其国内价格更接近"封闭系统均衡",适合拟合马来西亚这种相对封闭的市场;出口占比高(贸易型)的国家,价格被国际市场拉着走,反而不适合。脚本顶部 target_country 列了 20 个候选国(澳大利亚、西班牙、摩洛哥、英国、波兰、法国……),country selection.py 的输出就是从这个池子里挑出最合适的控制组。

脚本同时还算了"出口金额占比"(用价格加权),与"出口量占比"交叉验证:

#country selection.py:217
value['Value'] = np.multiply(value['Value_x'], value['Value_y'])  # 产量×单价=总价值
​

两个口径都算一遍,稳健性更高。这是政策评估文献里的标准做法。

算法与代码精读

prepare:五表 join 成大宽表

cleansing 的核心函数 prepare(target_land, target_prod, target_prix)(L19-80)做的事是:按年份循环,把生产、面积、价格、单产逆(yield_i = 面积/产量)四张子表两两 merge,拼成一张大宽表 D[year]。

#cleanse data.py:54-66(节选)
for currentyear in range(beginyear, endyear):
    temp1 = target_prod[...].merge(target_area[...], on=['Item','Year'], how='outer')
    temp2 = target_prix[...].merge(target_yield_inverse[...], on=['Item','Year'], how='outer')
    # 重命名列:Value_x→production, Value_y→area, Value_x→price, Value_y→yield_i
    temp1.columns = temp1.columns.str.replace('Value_x','production')
    temp1.columns = temp1.columns.str.replace('Value_y','area')
    ...
    data = temp1.merge(temp2, on=['Item','Year','type','lifespan'], how='outer')
    D[currentyear] = data
​

注意 yield_i(yield inverse)——它不是单产,而是单产的倒数 = 每吨产量需要的面积。这是为了第 01 节讲的等式约束 Σ Q_i × yield_i = 总可耕地,用倒数才能把"产量"换算成"占地"。

多年生寿命约束的计算

prepare 函数后半段(L68-78)给每年每作物算 eco lifespan:

#cleanse data.py:69-78
for currentyear in range(beginyear, endyear):
    eco_lifespan = [default_discount for _ in range(len(D[currentyear]))]   # 默认 0.8
    indices = D[currentyear]['lifespan'].dropna().index
    perennial = D[currentyear]['lifespan'].dropna().apply(
        lambda x: (x*eco_coeff - 1) / (x*eco_coeff)                         # (ωL-1)/(ωL)
    )
    for i in indices:
        eco_lifespan[i] = round(perennial.loc[i], 4)
    D[currentyear]['eco lifespan'] = eco_lifespan
​

逻辑分三步:① 默认所有作物都是一年生,用 default_discount=0.8(惯性系数 ψ);② 对 mapping.csv 里有寿命的多年生作物,改用 (ωL−1)/(ωL) 算下界系数;③ round 到 4 位小数写回。

注意公式 (ωL−1)/(ωL) 算出来是负数(因为 ωL 通常远大于 1,比如 ωL=21 时结果是 −20/21 ≈ −0.95)。这正好契合 forecast.py 里 inequality_value=cvxopt.matrix(np.multiply(np.multiply(area_hist,eco_lifespan),-1)) 的写法——eco_lifespan 已经是负数,再乘 −1 就成了正的面积阈值。这是作者为了把 ≥ 约束塞进 CVXOPT 的 ≤ 标准形式留下的符号痕迹,读代码时要小心。

合成控制价格的填充

prepare 返回的 grand 表里,6 种作物的 price 列还是 NaN。L300-308 做填充:

#cleanse data.py:301-308
fillnull = dict(zip(
    grand['Item'][grand['price'].isnull()].unique(),
    [malay_tobacco, malay_ginger, malay_lettuce,
     malay_orange, malay_maize, malay_sugarcane]
))
for i in fillnull:
    indices = grand[grand['Item']==i][grand['price'].isnull()].index.tolist()
    for j in indices:
        grand.at[j, 'price'] = fillnull[i][indices.index(j)]
​

按作物名 → 合成价格列表的字典,逐 NaN 填充。填完后导出 grand.csv,供下一节 estimate demand.py 用。

两个工程坑

坑①:global 变量满天飞。整个脚本顶层的 beginyear / endyear / eco_coeff / default_discount 都是 global,函数体内还出现 global D 声明(prepare 函数 L50)。这是 Jupyter Notebook 风格的遗留——cell 之间的状态全靠 global 传递。改造成模块化代码时要把这些 global 显式当参数传进去。

坑②:旧版 pandas append。脚本里多次用 malay_gdp.append(gdp_extra)(L517)、malay_land.append(land_extra)(L548)、malay_prix.append(...).append(...)(L202)。pandas 1.4 起 DataFrame.append 已废弃,新版会报 AttributeError: 'DataFrame' object has no attribute 'append'。改法是换成 pd.concat([malay_gdp, gdp_extra])。作者当年(2020)写的代码在 pandas 1.x 跑得欢,现在直接跑会全红。

🎯 反直觉发现:合成控制法原本是评估"政策效果"的工具,Smart Farmers 把它拿来"造数据"。在 Abadie 原文里,合成巴斯克 vs 真实巴斯克的 GDP 差,衡量的是恐怖活动的经济破坏;在 Smart Farmers 里,合成马来西亚 vs(未来真实的)马来西亚价格差,衡量的是模型预测误差。同一个数学(加权 OLS 拟合 + 外推预测),用途从"因果推断"变成了"数据插补"。这是跨界借用的典范——好方法的灵魂是数学结构,不是它最初被发明出来的那个语境。

关键技巧

1. join key 三件套:Area(国家)、Item(作物)、Year(年份)是 FAO 五表的通用 key,所有 merge 都靠它。
2. subtotal 必须剔除:"Cereals, Total"这种汇总类会和细分类重复计算,check consistency.py 的 exclude 清单就是干这个的。
3. yield 的倒数 yield_i:建模时用"每吨所需面积"而不是"每亩产量",是为了让等式约束 Σ Q × yield_i = 总地 形式干净。
4. SCM 的角色翻转:任何"加权拟合 + 外推"的方法都可以从"评估处理效应"变成"补缺失值",前提是处理组和控制组在"非处理期"行为相似。
5. Google 首页取均值:粗糙但诚实的取数规则。学术上不优雅,工程上够用,关键是显式声明让读者可复现可质疑。
6. 旧 pandas API 迁移:看到 .append( 一律改 pd.concat([...]),看到 .append(...).append(...) 改 pd.concat([a, b, c])。

本节要点

1. FAO 数据库 14GB,五个核心表(Production / Prices / LandUse / Population / Macro),用 Area+Item+Year 三键 join。
2. 马来西亚 61 类作物 → 剔 subtotal + 剔 <1% 用地 → 30 类,时间窗 2013-2018 回测 + 2019-2025 前向。
3. 多年生作物寿命表 mapping.csv 用"Google 首页取均值",下界系数 (ωL−1)/(ωL)(ω=0.7,L=作物寿命)。
4. 合成控制法补 6 种缺价格作物(姜/生菜/玉米/橙/甘蔗/烟草),灵感来自 Abadie 2003 巴斯克研究,角色从"评估政策"翻转为"补缺失值"。
5. country selection.py 用"出口量/产量"和"出口金额/产值"两个口径筛控制组国家。
6. 工程坑:global 变量满天飞(Jupyter 遗留)、malay_gdp.append 是 pandas 旧 API,新版要换 pd.concat。

下一节,我们进入 Smart Farmers 的计量经济学心脏——estimate demand.py 怎么用工具变量(2SLS)处理"需求不可观测"的难题,把定价系数 α、人口系数 β、GDP 系数 γ 估计出来,塞进 forecast.py 的 QP。