第 7 章 · 03 需求估计与工具变量 2SLS

第 7 章 · 03 需求估计与工具变量 2SLS

本节摘要:本节精读 estimate demand.py(约 326 行),Smart Farmers 的计量经济学心脏。上一节我们凑齐了 30 类作物的价格、产量、人口、人均 GDP,但第 01 节的定价机制 P=P_eq+α(D−Q)+ε 里还藏着一个幽灵——国内需求 D 不可观测。你拿不到"马来西亚每年消费了多少吨甘蓝"这种数据。本节讲作者怎么用**工具变量(Instrumental Variable)+ 两阶段最小二乘(2SLS)**绕过这个难题:人口与人均 GDP 对作物价格没有直接因果,只能经由需求间接影响,是完美的工具变量。具体到代码,作者用 constrained_ols 函数把带非负约束的最小二乘写成 QP(又是 CVXOPT!),串联两个回归(需求模型 + 定价模型),最终估计出每个作物的三个系数 α(定价)/ β(人口)/ γ(GDP),输出 tres_grand.csv 供下一节 forecast.py 用。本节是全书计量味最浓的一节,但我们会把它讲得直觉优先。

涉及源码:原项目 Smart Farmers project/estimate demand.py(约 326 行),输出 tres_grand.csv(下一节 forecast.py 的主输入)。

⚠️ 注意:本节的数学有两个层次——一是工具变量为什么能解决"需求不可观测",二是 constrained_ols 为什么用 QP 而不是普通 OLS。两个层次都搞懂了,你才算真正读懂 Smart Farmers。如果你没学过计量经济学,本节会从直觉讲起,不必先去看教材。

学习目标

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

1. 说清"需求 D 不可观测"为什么是定价机制的死穴。
2. 解释工具变量的两个充要条件(相关性 + 外生性),以及人口 + GDP 为何同时满足。
3. 复述 2SLS 的两阶段流程(第一阶段拟合 D,第二阶段把拟合的 D 代入定价模型)。
4. 读懂 constrained_ols 怎么把带非负约束的 OLS 写成 QP(P=XᵀX, q=−Xᵀy)。
5. 解释为什么常数项不加非负约束(inequality_coeff[0,0]=0)。
6. 说清 tres_grand.csv 在整个 pipeline 里的位置(承接 cleansing,下接 forecast)。

金融/数学原理

幽灵变量:不可观测的国内需求

回到第 01 节的定价机制:

P_i = P_eq,i + α_i × ( D_i − Q_i ) + ε_i
​

这里 Q_i(产量)我们有数据——FAO Production 表;P_i(价格)我们有数据——FAO Prices 表,缺的用合成控制法补;但 D_i(国内需求)呢?没有任何数据库直接给"马来西亚每年消费了多少吨甘蓝"。这就是定价机制里的幽灵变量。

如果不处理这个幽灵,定价模型就退化不成样子。你可能会想用一些代理变量——比如把"产量 − 出口量"当需求。但 FAO 的出口数据有 4GB 之巨,而且很多作物的出口口径混乱(有的按鲜重,有的按干重),这条路在工程上很痛苦。

作者选择了另一条路——用工具变量绕过需求的直接观测。

工具变量的直觉:绕路,不硬刚

工具变量(Instrumental Variable,IV)是计量经济学的镇山之宝,专门处理"自变量里掺了不可观测误差"的情形。它的核心是换一条因果路径绕过去。

考虑一个简化版定价模型:

P = α·D + 扰动
D = 真实需求 + 不可观测误差
​

如果我们直接拿 D 回归 P,不可观测误差会污染估计。工具变量的思路是:找一个变量 Z,它满足两个条件:

满足这两个条件的 Z 就是工具变量。找到 Z 后,用 两阶段最小二乘(2SLS):

• 第一阶段:用 Z 回归 D,得到 D 的拟合值 D̂(把 D 里被不可观测误差污染的部分"洗掉");
• 第二阶段:用 D̂ 回归 P,得到的系数就是干净的因果效应。

人口 + GDP:完美的工具变量对

作者在 README 里点出:

Both population and GDP per capita make brilliant instrumental variables. In layman's terms, population and GDP per capita have no direct power over the crop price. They can only influence the crop price indirectly through the change of demand.

直觉非常清楚:

• 人口↑ → 吃饭的嘴变多 → 需求 D↑ → 价格 P↑:人口对价格只有这一条路径,相关性成立;
• 人均 GDP↑ → 大家更有钱 → 需求 D↑(尤其是高价值作物)→ 价格 P↑:GDP 对价格也只有这一条路径,外生性成立;
• 反过来,人口或 GDP 直接决定甘蓝价格?不可能。人口再多,如果大家不吃甘蓝,价格也不动。这就是外生性。

所以人口和 GDP 是天然的工具变量。Smart Farmers 里 β 对应人口系数、γ 对应 GDP 系数、α 是定价系数,这三个系数正是下一节 forecast.py 的 D[currentyear]['alpha'/'beta'/'gamma']。

Angrist & Krueger (2001) 的传承

README 文献第 2 篇是 Angrist & Krueger (2001) 的工具变量综述。这篇综述总结了用工具变量分析自然实验和随机实验数据的经典案例,核心论点是:工具变量能在不可观测变量存在的情况下识别因果。Smart Farmers 对此的应用很纯粹——需求 D 就是那个不可观测变量,人口 + GDP 是工具,2SLS 是识别手段。

💡 核心心法:工具变量的价值不是"让回归更准",而是"在没有随机对照实验的情况下,从观测数据里抢救出因果"。这是观察性研究里最接近"准实验"的工具。Smart Farmers 没法做"随机给马来西亚人发甘蓝"的实验,但人口和 GDP 的自然变异就是天然的实验条件。

为什么用 QP 而不是普通 OLS:非负约束

这是本节的第二个数学层次。普通最小二乘(OLS)的解是闭式的:β̂ = (XᵀX)⁻¹ Xᵀy。但作者偏要用 CVXOPT 解 QP 来做回归,为什么?

因为这里需要带非负约束——经济学上,价格对供给的敏感度 α 应该是正数(供给过剩价格跌,供给不足价格涨,α>0);人口越多需求越大,β>0;GDP 越高高价值作物需求越大,γ>0。无约束 OLS 可能估出负的 α/β/γ,这违反经济学直觉。

带非负约束的 OLS 是一个标准 QP:

min   ||Xβ − y||²    s.t.   β ≥ 0
​

展开目标函数:||Xβ − y||² = βᵀXᵀXβ − 2yᵀXβ + yᵀy,常数项 yᵀy 不影响最优解,所以等价于:

min   (1/2)βᵀ(2XᵀX)β + (−2Xᵀy)ᵀβ    s.t.   −β ≤ 0
​

这正是 CVXOPT 的 QP 标准形式 min (1/2)xᵀPx + qᵀx s.t. Gx≤h,其中 P=2XᵀX(可省略因子 2,因 CVXOPT 接受任意正定 P)、q=−Xᵀy、G=−I、h=0。

算法与代码精读

create_xy:构造回归设计矩阵

estimate demand.py 的 create_xy(target_crop, grande, malay_gdp, malay_pop) 把价格 y 和自变量 X 拼好:

#estimate demand.py:22-35
def create_xy(target_crop, grande, malay_gdp, malay_pop):
    y = grande['price'][target_crop]                                # 价格(因变量)
    x = pd.concat([malay_gdp['Value'], malay_pop['Value'],
                   grande['production'][target_crop]], axis=1)      # GDP, 人口, 产量
    x = sm.add_constant(x)                                          # 加常数项
    # 把产量列取负号(因为定价机制里是 D−Q,Q 前是负号)
    x[target_crop.lower()] = x[target_crop].apply(lambda x: -1*x)
    del x[target_crop]
    return x, y
​

这里有个关键的符号 trick:production 列被取了负号,然后列名改成小写。为什么?回到定价机制 P = P_eq + α(D − Q),展开 P = P_eq + α·D − α·Q。回归设计矩阵里 D 被人口和 GDP 替代(工具变量),Q 前面带负号,所以 X 里产量列乘 −1。这样回归得到的 α 直接是正数(经济学正确)。

注意 sm.add_constant 加的常数项对应 P_eq(均衡价格基准)。

constrained_ols:QP 解带约束的最小二乘

本节的核心函数 constrained_ols(x, y)(L76-97):

#estimate demand.py:76-97
def constrained_ols(x, y):
    linear_coeff    = cvxopt.matrix(-1*np.mat(y.tolist())*np.mat(x)).T   # q = -X^T y
    quadratic_coeff = cvxopt.matrix(np.mat(x).T*np.mat(x))               # P = X^T X

    # 不等式约束:β ≥ 0 转成 -β ≤ 0
    inequality_coeff = cvxopt.matrix(0.0, (len(x.columns), len(x.columns)))
    inequality_coeff[::len(x.columns)+1] = -1                            # 对角线填 -1(即 -I)
    # 常数项不加约束!
    inequality_coeff[0, 0] = 0
    inequality_value = cvxopt.matrix([0.0 for _ in range(len(x.columns))])

    cvxopt.solvers.options['show_progress'] = False
    ans = cvxopt.solvers.qp(P=quadratic_coeff, q=linear_coeff,
                            G=inequality_coeff, h=inequality_value)['x']
    return ans
​

逐行拆:

1. quadratic_coeff = np.mat(x).T * np.mat(x):就是 XᵀX。CVXOPT 要求 P 正定(或半正定),XᵀX 天然对称半正定,满足。
2. linear_coeff = -1 * np.mat(y) * np.mat(x):这是 −yᵀX,转置后就是 q = −Xᵀy。注意代码里写法是 y·x,然后 .T 转置,等价。
3. inequality_coeff[::len+1] = -1:这是用一维索引填对角阵的 CVXOPT 惯用法。len+1 是 n×n 矩阵里对角元素的步长(0, n+1, 2n+2, ...),填 −1 就成了 −I 单位阵。
4. inequality_coeff[0,0] = 0:**常数项不加约束!**因为常数项对应 P_eq(均衡价格),价格本来就该是正数,加非负约束是冗余;更重要的是,如果常数项也被压到 ≥0,会限制拟合的灵活性。这是作者的细心。
5. cvxopt.solvers.qp:标准 QP 求解,返回 α/β/γ/常数四元组(顺序按 X 的列:const, GDP, pop, production_neg)。

get_params:对所有作物跑一遍 constrained_ols

#estimate demand.py:104-136(节选)
def get_params(crops, grande, malay_gdp, malay_pop, viz=False):
    D = {}
    for target_crop in crops:
        x, y = create_xy(target_crop, grande, malay_gdp, malay_pop)
        ans = constrained_ols(x, y)
        if viz:
            forecast = np.mat(ans).T * np.mat(x).T
            forecast = forecast.ravel().tolist()[0]
            # 画 Est vs Act 对比图
            ...
        D[target_crop] = list(ans)
    return D
​

对 30 类作物逐个跑 constrained_ols,把每类的系数四元组存进字典。viz=True 时画预测价 vs 真实价对比图(README 里的 cabbage/cocoa/coconut/mango/rubber/oil palm regression 图就是这么生成的)。

输出:列重命名 → 合并到 grand → tres_grand.csv

L217-237 把字典转成 DataFrame,重命名列,合并回 grand:

#estimate demand.py:217-237
output = pd.DataFrame(columns=D2.keys())
for i in D2:
    output[i] = D2[i]
output = output.T
output.columns = ['constant', 'gamma', 'beta', 'alpha']   # 注意顺序!
output.index.name = 'Item'
output.reset_index(inplace=True)

# 合并回 grand
tres_grand = grand.merge(output, on='Item', how='left')
tres_grand.to_csv('tres_grand.csv', index=False)
​

列顺序的坑:output.columns 按字段命名为 [constant, gamma, beta, alpha]。回忆 create_xy 里 X 的列顺序是 [const, GDP, pop, production_neg],所以 gamma 对应 GDP、beta 对应人口、alpha 对应产量(定价系数)。命名时要小心别搞反 β 和 γ。

tres_grand.csv 就是下一节 forecast.py 的主输入。文件名 tres_grand(法语"非常大的 grand")暗示它在 grand 的基础上又加了三列回归系数,变成了一个"更大的 grand"。

两个回归的串联:需求模型 + 定价模型

虽然 estimate demand.py 的代码只显式跑了一个回归(定价模型),但 README 的理论框架里其实是两个回归的串联——这正是 2SLS 的两阶段:

阶段一(需求模型):   D = β·人口 + γ·GDP + ε_D
阶段二(定价模型):   P = P_eq + α·(D̂ − Q) + ε_P
​

代码里的简化是:直接把人口和 GDP 当作 D 的代理塞进定价回归,跳过了显式估 D̂ 这一步。这在数学上等价于 2SLS(因为线性模型的 2SLS 可以证明等价于"把工具变量直接当自变量,但只解释被工具化的部分")。constrained_ols 估出的 β/γ 既是需求模型的系数,也是定价模型里 D 的分量系数——一箭双雕。

上图是 README 里甘蓝(cabbage)的回归拟合结果。紫色(估)捕捉了黄色(真)的整体趋势但波动更小——这是约束回归的典型特征:非负约束压住了过拟合,拟合线更平滑。README 评价:

In the actual regression stage, the in-sample data is a nice fit. The purple line captures the overall trend but has a smaller volatility than the yellow line.

实证结果速览(6 种作物)

README 的 Empirical Result 段展示了 6 种作物的回测图,每种代表一种挑战:

整体准确率约 65%(README 原话"a useful tool on around 65% of the crops")。作者宽慰:"In quantitative trading, we can churn out a Gulfstream G6 from a factor with odds at 55%"——量化交易里 55% 胜率的因子就能买湾流 G6,65% 已经很好了。

🎯 反直觉发现:constrained_ols 用 QP 解 OLS,看起来是杀鸡用牛刀,实则是为了"用约束守卫经济学含义"。普通 OLS 是闭式解 β̂=(XᵀX)⁻¹Xᵀy,一行代码搞定;这里偏要绕道 CVXOPT 解 QP,工程上慢得多。但代价换来的是 α/β/γ 的非负性——无约束 OLS 在小样本(马来西亚 6 年数据)下很容易估出负的价格系数,这在经济学上荒谬。约束回归的本质不是"拟合更好",而是"保证估出来的系数有意义"。这是量化建模里反复出现的主题:可解释性 > 拟合优度。

关键技巧

1. 设计矩阵的符号 trick:定价机制里 Q 前是负号,设计矩阵里产量列直接乘 −1,这样回归系数天然为正。比"先估正号再翻号"更不容易出错。
2. CVXOPT 填对角阵的惯用法:matrix[::n+1] = -1 用一维索引填 n×n 对角阵,比 np.diag([-1]*n) 更紧凑(但可读性差,新手慎用)。
3. 常数项不加约束:inequality_coeff[0,0]=0 是点睛之笔。常数项对应 P_eq,加非负约束既冗余又损失灵活性。
4. 2SLS 的代码等价:线性模型下,2SLS 等价于"直接把工具变量当自变量回归"——所以代码里只跑一个回归,数学上却完成了两阶段。
5. viz 标志位:get_params 的 viz=False 默认关画图,跑全样本时不开 viz 避免刷屏;调试单个作物时 viz=True 看拟合。
6. 文件命名 tres_grand:法语"非常大的 grand",暗示它是 grand 表加了 α/β/γ 三列后的"升级版"。作者的命名幽默感。

本节要点

1. 定价机制 P=P_eq+α(D−Q)+ε 的死穴是需求 D 不可观测——没有任何数据库直接给"年消费量"。
2. 工具变量 IV + 2SLS 是绕过不可观测变量的标准手段。人口和 GDP 对价格无直接因果(只经需求间接影响),是完美的工具变量对。
3. constrained_ols 把带非负约束的 OLS 写成 QP:P=XᵀX, q=−Xᵀy, G=−I, h=0,常数项不加约束。
4. 估计出每个作物的 α(定价系数)、β(人口系数)、γ(GDP 系数),输出 tres_grand.csv。
5. 6 种作物回测整体准确率约 65%,甘蓝最贴合(纯内销),芒果最糟(进口替代),油棕产量预测误差最小但价格最难投机。
6. 2SLS 的价值不是更准,而是在不可观测变量下识别因果;constrained_ols 的价值不是更好拟合,而是保证系数经济学含义正确。

下一节,我们终于要解那个压轴的 QP 了——forecast.py 怎么把 α/β/γ 组装成 quadratic_coeff 和 linear_coeff,用 cvxopt.solvers.qp 一次解出 30 维最优产量,以及用 Nelder-Mead 找最优成本初值。还有全书最深的教训:前向测试为什么打脸回测。