第 7 章 · 04 凸优化求解与前向测试

第 7 章 · 04 凸优化求解与前向测试

本节摘要:本节精读 forecast.py(约 618 行)和 check consistency.py(约 149 行),Smart Farmers 的工程落点与思想高潮。前三节搭好了所有零件——本节把它们组装成一台会跑的机器。forecast.py 的 get_production 把上一节估出的回归系数 α/β/γ 组装成 CVXOPT 的 quadratic_coeff(α 对角阵)和 linear_coeff(由历史价格、历史产量、Δ人口、ΔGDP、成本初值拼接),用 cvxopt.solvers.qp 解出 30 维最优产量;Nelder-Mead(scipy.optimize.minimize)在 costfunction 上随机搜索最优成本初值,cost_optimal 那 30 个硬编码数字就是搜索结果。前向测试外接 IMF 人均 GDP(capita.csv)+ UN 人口预测,假设农地年增 1%/5%,对 2019-2025 情景预测,再用 Bursa Malaysia 棕榈油 CPO 期货做前向验证。但本节最深刻的不是求解,而是前向测试打脸回测——模型预测 2018-2020 backwardation、2021+ contango,实际却被 COVID 的 W 型复苏砸碎。引 MSCI 研究主管 Dimitris Melas 那句"I have never seen a bad backtest"。最后讲 Ricardian 推广版:把 Pareto Optimal 升级成 Nash Equilibrium,禁止国家间竞争合谋(点名 2020 俄沙石油战),作者务实承认"收集全部数据要十年",选择不修未坏之物。

涉及源码:原项目 Smart Farmers project/forecast.py(约 618 行),Smart Farmers project/check consistency.py(约 149 行),外接 capita.csv(IMF 人均 GDP 预测)。

⚠️ 注意:本节代码量最大,涉及 CVXOPT 的 QP 标准形式转换、scipy 的 Nelder-Mead、双轴可视化、前向测试的多个外生假设。建议先把第 01 节的目标函数推导和第 03 节的 α/β/γ 含义复习一遍,本节的代码才看得懂。最深刻的反直觉发现放在"前向测试"段,那是全章的思想高潮。

学习目标

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

1. 逐行讲清 get_production 怎么把 α/β/γ 组装成 CVXOPT 的 P/q/G/h/A/b。
2. 解释轮作上下界怎么用 upperblock/lowerblock 块矩阵拼成一个 G。
3. 说清 Nelder-Mead 在本项目的角色——找最优成本初值,不是直接找最优产量。
4. 复述 costfunction 的 SSE 归一化形式 Σ((actual−est)/actual)² 及其动机。
5. 讲清前向测试的三个外生输入(IMF GDP / UN 人口 / 农地年增假设)。
6. 解释为什么前向测试远不如回测——"上帝视角 vs 现实"、"COVID W 型复苏"、"I have never seen a bad backtest"。
7. 复述 Ricardian 推广版的两条加强假设(Pareto Optimal → Nash Equilibrium,数据可用性)。

金融/数学原理

回顾:我们已经有了什么

经过前三节,我们手上已经有:

• 30 类作物的价格、产量、面积、单产逆(yield_i)、寿命(多年生)、eco lifespan(下界系数);
• 每类作物的回归系数 α(定价)、β(人口)、γ(GDP),存在 tres_grand.csv;
• 马来西亚的人口(malay_pop.csv)、人均 GDP(malay_gdp.csv)、可耕地(malay_land.csv)时间序列;
• 一个清晰的目标函数:最大化利润 = Σ(P−C)·Q,其中 P=P_eq+α(D−Q),所以利润对 Q 是凹二次型。

本节的任务:给定上述输入,求出每年每类作物的最优产量 Q*,然后把这套机制外推到 2019-2025。

CVXOPT 的 QP 标准形式

CVXOPT 要求所有 QP 都塞进这个形式:

minimize   (1/2) xᵀ P x + qᵀ x
subject to   G x ≤ h
             A x = b
​

把 Smart Farmers 的模型对号入座:

forecast.py 的 get_production 就是按这张表逐行填矩阵。

三个外生假设:前向测试的命门

回测时,所有输入都来自历史数据(2013-2018)——人口、GDP、可耕地、单产全是真实的,我们站在 2020 年回看,拥有"上帝视角"。但前向测试(2019-2025)做不到,因为这些输入的"未来值"必须预测。作者用了三个假设:

这三个假设是前向测试的命门——任何一个偏离现实,预测就崩。后文会看到,COVID 把 GDP 和人口预测双双打废,这是失败的主因。

算法与代码精读

get_production(二次版):核心 QP 组装

forecast.py 同时定义了线性版(L55-86)和二次版(L93-150),二次版覆盖线性版,是实际使用的。逐段拆二次版:

#forecast.py:93-123(节选)
def get_production(initial_guess):
    # ① 取历史值
    price_hist       = np.mat(D[currentyear-1]['price']).T          # P_eq 近似(去年价)
    alpha            = np.diag(D[currentyear]['alpha'])             # 定价系数对角阵
    production_hist  = np.mat(D[currentyear-1]['production']).T     # 去年产量 Q_{t-1}
    delta_pop        = malay_pop[...currentyear...item() - ...currentyear-1...item()]   # Δ人口
    beta             = np.mat(D[currentyear]['beta']).T             # 人口系数
    delta_gdp        = malay_gdp[...currentyear...item() - ...currentyear-1...item()]   # ΔGDP
    gamma            = np.mat(D[currentyear]['gamma']).T            # GDP 系数

    # ② 线性项 q
    linear_coeff = (-price_hist
                    + (alpha*production_hist*-1)       # 注意这里是 -α·Q_{t-1}
                    - delta_pop*beta                    # -Δ人口·β
                    - delta_gdp*gamma                   # -ΔGDP·γ
                    + np.mat(initial_guess).T)          # +成本初值(对应 -C·Q 里的 C)
    linear_coeff = cvxopt.matrix(linear_coeff)

    # ③ 二次项 P
    quadratic_coeff = cvxopt.matrix(alpha)              # 就是 α 对角阵
​

关键观察:quadratic_coeff = cvxopt.matrix(alpha)——目标函数的二次项矩阵就是 α 对角阵。这正是第 01 节推导的"价格内生 → 二次型"在代码里的兑现。

linear_coeff 的符号有点绕,需要回到目标函数展开:

利润 = Σ [ P_eq + α(D − Q_{t}) − C ] × Q_{t}
     = Σ [ P_eq + α·D − α·Q_{t} − C ] × Q_{t}
     = Σ [ (P_eq + α·D − C)·Q_{t} − α·Q_{t}² ]
​

最大化利润等价于最小化 −利润,所以 CVXOPT 的 q 是 −利润的线性部分:

q = −(P_eq + α·D − C)
  = −P_eq − α·D + C
​

其中 α·D 这一项要用工具变量展开 α·D = α(β·Δpop + γ·Δgdp + 历史 D)≈ β·Δpop + γ·Δgdp + α·Q_{t-1}(用历史产量近似历史需求)。所以:

q = −P_eq + α·Q_{t-1} − Δpop·β − ΔGDP·γ + C(initial_guess)
​

完全对应代码里 linear_coeff 的五项。注意 α·Q_{t-1} 在代码里写成 alpha*production_hist*-1(因为整体符号取反,代码里 production_hist 已经是 +号,再乘 −1 就成了 −α·Q_{t-1},与公式中的 +α·Q_{t-1} 在最终 −利润 形式下一致)。

约束矩阵:G/h 块拼接

#forecast.py:125-143(节选)
    # 等式约束:总可耕地
    equality_coeff = cvxopt.matrix(D[currentyear]['yield_i']).T      # A = yield_i 向量
    equality_value = cvxopt.matrix(malay_land['Value'][...].item()*1000)   # b = 总可耕地(×1000 单位换算)

    # 不等式约束:轮作上下界
    area_hist      = D[currentyear-1]['area']                        # 去年面积
    eco_lifespan   = D[currentyear]['eco lifespan'].tolist()         # 多年生下界系数

    upperblock = np.diag(D[currentyear]['yield_i'])                  # 上界对角块(正)
    lowerblock = np.diag(D[currentyear]['yield_i'].apply(lambda x: -1*x))   # 下界对角块(负)

    # 垂直拼接成 60×30 的 G
    inequality_coeff = cvxopt.matrix(np.append(upperblock, lowerblock, axis=0))
    inequality_value = cvxopt.matrix(
        (area_hist*upperbound).append(np.multiply(np.multiply(area_hist,eco_lifespan),-1))
    )
​

三个细节:

1. 等式约束 A x = b:A 是单产逆向量(yield_i),b 是总可耕地面积。*1000 是单位换算(FAO 单位与 malay_land 单位差 1000)。
2. 上下界块拼接:upperblock 是正对角阵(实现 ≤ 上界),lowerblock 是负对角阵(实现 ≥ 下界,通过 −x ≤ −阈值)。两个块垂直拼成 60 行的 G。
3. inequality_value 也拼:(area_hist*upperbound).append(...) 把上界阈值和下界阈值拼成一个 60 维的 h。upperbound=1.2(轮作上界,允许面积年增 20%),下界用 area_hist × eco_lifespan × −1(eco_lifespan 对一年生是 −0.8,多年生是 (ωL−1)/(ωL),都是负数,再乘 −1 成正阈值)。

注意 (area_hist*upperbound).append(...) 又是旧 pandas API,新版要换 pd.concat。

get_ans:CVXOPT qp 的薄包装

#forecast.py:39-48
def get_ans(quadratic_coeff, linear_coeff, inequality_coeff,
            inequality_value, equality_coeff, equality_value):
    cvxopt.solvers.options['show_progress'] = False                  # 关闭迭代日志
    ans = cvxopt.solvers.qp(P=quadratic_coeff, q=linear_coeff,
                            G=inequality_coeff, h=inequality_value,
                            A=equality_coeff, b=equality_value)['x']
    return list(ans)
​

一行 cvxopt.solvers.qp 解出 30 维最优产量。CVXOPT 内部用内点法,30 维 QP 解起来毫秒级。['x'] 取解向量,转 list 返回。

costfunction + Nelder-Mead:找最优成本初值

这是 forecast.py 最巧妙也最容易被忽略的一段。问题:linear_coeff 里有一项 initial_guess(对应成本 C),这个 C 怎么定?

作者的方案是用 Nelder-Mead 在成本空间里搜索,让模型预测的产量与真实产量误差最小:

#forecast.py:157-166
def costfunction(initial_guess):
    estimate = get_production(initial_guess)                         # 用当前成本猜解 QP
    actual   = D[currentyear]['production'].tolist()                 # 真实产量
    cost = sum(np.power(np.divide(np.subtract(actual, estimate), actual), 2))   # 归一化 SSE
    return cost
​

costfunction 是一个"QP 嵌套"——外层成本函数调用内层 QP 求解。归一化形式 Σ((actual−est)/actual)² 用实际值除一下,避免大量级作物(油棕几千万吨)淹没小量级作物(烟草几千吨)的贡献。

#forecast.py:173-199(节选)
def ls_estimate(initial_guess, diagnosis=True):
    if diagnosis:
        sse_original = costfunction(initial_guess)
        print(f'Initial SSE: {round(sse_original,2)}')
    # Nelder-Mead 搜索最优成本
    least_square = scipy.optimize.minimize(costfunction,
                                           x0=(initial_guess),
                                           method='Nelder-Mead')
    if least_square.success:
        if diagnosis:
            sse = costfunction(least_square.x)
            print(f'Result SSE: {round(sse,2)}')
        return least_square.x
    else:
        print(least_square)
        return [0]
​

Nelder-Mead 是无梯度优化算法(单纯形法),适合这种"目标函数不可导、调用一次很贵"的情形。每次迭代都要解一次 QP,所以这个搜索计算很重——这就是 find_init 要做"随机搜索 + 多次试"的原因。

#forecast.py:206-232(节选)
def find_init(num=10):
    global currentyear
    dic = {}
    for _ in range(num):
        # 随机生成成本初值(用价格乘 [0,1] 随机数)
        initial_guess = pd.Series([i*rd.random() for i in D[currentyear]['price']])
        ans = ls_estimate(initial_guess, diagnosis=False)
        if len(ans) > 1:
            sample = []
            # 在所有回测年份上算平均 SSE
            for currentyear in range(beginyear, endyear):
                cost = costfunction(ans)
                sample.append(cost)
            dic[np.mean(sample)] = initial_guess
    return dict(sorted(dic.items()))
​

find_init 跑 num=10 次随机初值,每次让 Nelder-Mead 优化出对应的最优成本,然后在所有回测年份上算平均 SSE,把"SSE → 对应初值"排序返回。SSE 最小的那个初值就是全局最优成本猜测。

cost_optimal:30 个硬编码数字

forecast.py:361-390 有一段突兀的 30 个硬编码数字:

#forecast.py:361-390
cost_optimal = [46.45261724887184,
 25.771497795481803,
 8.04132847940984,
 ...
 27.02920135710488]
​

这就是作者用 find_init 跑出来、再手动固化下来的最优成本初值。30 维对应 30 类作物的单位成本(美元/吨)。这种"硬编码搜索结果"的做法不优雅(理想是脚本自动跑),但工程上很务实——Nelder-Mead 搜索很慢,固化结果后回测和前向测试都能直接用。

回测主循环:6 年预测

#forecast.py:430-467(节选)
X = {}
for currentyear in range(beginyear, endyear):        # beginyear=2013, endyear=2019
    x = get_production(cost_optimal)                  # 用 cost_optimal 解 QP
    X[currentyear] = x

# 画每作物 Est vs Act
for ii in range(len(X[currentyear])):
    Y = [D[i]['production'][ii] for i in range(beginyear, endyear)]
    fig = plt.figure(figsize=(10,5))
    ...
    plt.plot(range(beginyear, endyear), [X[i][ii] for i in X], label='Est', color='#ef6466')
    plt.plot(range(beginyear, endyear), Y, label='Act', color='#9f65a2')
    plt.title(D[currentyear]['Item'][ii]+' Production')
​

6 年循环,每年解一次 QP,画每作物的 Est(红)vs Act(紫)对比图。compute_price 函数(L257-284)用回归系数算预测价格,同样画对比图。

前向测试:2019-2025

#forecast.py:503
beginyear = 2019; endyear = 2025

# ① 接 IMF GDP 预测
gdp_extra = malay_gdp.iloc[:(endyear-beginyear)].copy()
gdp_extra['Value'] = capita['Mid Price'][str(beginyear):str(endyear)].tolist()   # IMF 预测
malay_gdp = malay_gdp.append(gdp_extra)

# ② 假设可耕地年增 1%
temp = [malay_land['Value'].iloc[-1]]
for i in range(len(land_extra)):
    temp.append(temp[-1]*1.01)                       # 年增 1%
land_extra['Value'] = temp
malay_land = malay_land.append(land_extra)

# ③ 前向预测循环
for currentyear in range(beginyear, endyear):
    D[currentyear] = D[currentyear-1].copy()
    D[currentyear]['production'] = get_production(cost_optimal)
    D[currentyear]['price']      = compute_price(np.mat(D[currentyear]['production']).T)
    D[currentyear]['area']       = np.multiply(D[currentyear]['production'], D[currentyear]['yield_i'])
​

注意每年 D[currentyear] = D[currentyear-1].copy()——前向时没有历史可参照,直接复用上一年的结构(单产逆、寿命、回归系数都假设不变),只有 production/price/area 每年重算。

前向验证:棕榈油 CPO 期货

estimate demand.py 后半段(L245-326)做棕榈油前向验证——把模型预测的 oil palm fruit FFB 价格与 Bursa Malaysia 的 CPO 期货.generic first 实际价格对比。

#estimate demand.py:277-282(节选)
# 用林吉特→美元换算(0.23)拼历史与期货预测,再取年均
temp = palm['Palm oil'][str(beginyear):].apply(lambda x: x*0.23).append(palm_futures['Palm oil'][:str(endyear+5)])
palmoil = temp.resample('1A').mean()
palmoil.index = [pd.to_datetime(str(i)[:5]+'01-01') for i in palmoil.index]
​

x*0.23 是林吉特→美元的近似汇率换算(CPO 期货以林吉特计价,FFB 价格以美元计价)。resample('1A').mean() 取年均,把日频期货转成年频对齐模型。

双轴可视化:oil palm vs palm oil

estimate demand.py:300-325 用 matplotlib 的 twinx() 画双轴图——左轴棕榈油期货价格,右轴油棕果模型价格(单位不同必须双轴):

fig = plt.figure(figsize=(10,5))
ax  = fig.add_subplot(111)
ax.set_xlabel('Date')
ax.set_ylabel('Palm Oil Generic 1st', color='#45ADA8')
ax.plot(palmoil.index, palmoil, color='#45ADA8', label='Palm Oil')
ax2 = ax.twinx()                                              # 双轴
ax2.set_ylabel('Oil Palm Fruit Price', color='#EC2049', rotation=270)
ax2.plot(palmoil.index[:len(oilpalm_act)], oilpalm_act,  color='#EC2049', label='Oil Palm Act')
ax2.plot(palmoil.index[len(oilpalm_est)-1:], oilpalm_est, color='#fea5c4', label='Oil Palm Est', linestyle='-.')
plt.title('Palm Oil vs Oil Palm')
​

check consistency.py:数据合理性体检

check consistency.py(149 行)是数据层的"体检脚本",做三件事:

1. subtotal 剔除的合理性:对比"按细分类求和"与"按 subtotal 类求和"的面积,验证 exclude 清单是否合理;
2. 筛选后面积 vs 全国可耕地:确保筛选出的 30 类作物面积之和不超过 malay_land 报告的 Cropland 总量;
3. 作物交集核验:a.intersection(b) 检查同时有"Area harvested"和"Production"数据的作物。

这是科研代码的好习惯——主 pipeline 之外另写一个 sanity check 脚本,不进生产但留作 review。

🎯 反直觉发现(本节思想高潮):前向测试远不如回测,模型预测被 COVID 的 W 型复苏打脸。回测时,作者站在 2020 年回看 2013-2018,所有输入(人口、GDP、可耕地)都是已发生的真实数据,等于"上帝视角"——所以回测漂亮,模型"看起来"在 65% 的作物上有效。但前向预测 2019-2025 时,作者用的是 IMF 和 UN 在 COVID 之前发布的人口与 GDP 预测,假设农地年增 1%、单产不变。模型给出"2018-2020 backwardation(近月 > 远月,现货偏紧)、2021+ contango(近月 < 远月,现货宽松)"的预测。现实是 COVID 在 2020 年砸了下来,各国 GDP 暴跌、人口增长受阻、供应链紊乱,棕榈油价格走出 W 型复苏——和模型的 backwardation/contango 节奏完全错位。作者引用 MSCI 研究主管 Dimitris Melas 的名言作为 Discussion 段的开篇:"I have never seen a bad backtest"。这句话道破量化研究的核心痛点——上帝视角的回测永远是好的,真正的考验是前向。如果 2018 年真的按模型做空棕榈油期货并持续滚动,2020 年基金会被赎回得精光。这个失败不是模型的错,是所有依赖外生预测的模型的共同命运——只要外生预测错,模型就错。Smart Farmers 的诚实在于把这次失败完整展示出来,而不是只秀回测。

关键技巧

1. QP 的 P 就是 α 对角阵:理论与代码的最紧对应。看到 quadratic_coeff = cvxopt.matrix(alpha) 这行,你应该立刻反应"这是价格内生带来的二次项"。
2. 块矩阵拼接实现多段约束:上界 + 下界拼成 60×30 的 G,是 CVXOPT 要求"单一不等式约束"的标准做法。
3. Nelder-Mead 找最优成本:外层无梯度优化 + 内层 QP,适合"目标函数不可导但可计算"的嵌套问题。每次迭代成本高,所以结果硬编码复用。
4. 归一化 SSE:Σ((actual−est)/actual)² 用实际值归一,避免量级悬殊的作物相互淹没。
5. 双轴可视化:不同单位的序列对比必须 twinx,否则大量级序列会把小量级压成一条平线。
6. check consistency 脚本:主 pipeline 之外另写 sanity check,是科研代码必备的"防呆"层。

本节要点

1. get_production 把 α/β/γ 组装成 CVXOPT 的 P(α 对角阵)/ q(5 项符号拼接)/ G(上下界块)/ h / A(yield_i)/ b(总可耕地)。
2. 约束三转换:最大化→最小化负目标、≥→≤(用 −1)、多段约束→块矩阵垂直拼接。
3. Nelder-Mead 在 costfunction 上搜索最优成本初值,cost_optimal 那 30 个硬编码数字是搜索结果。
4. costfunction 是 QP 嵌套——外层无梯度优化,内层每次解一次 QP,归一化 SSE 衡量产量预测误差。
5. 前向测试三外生假设:IMF GDP(capita.csv)、UN 人口、农地年增 1%、单产历史均值。
6. 棕榈油 CPO 前向验证:模型预测 2018-2020 backwardation、2021+ contango,实际被 COVID 的 W 型复苏打脸。
7. 思想高潮:"I have never seen a bad backtest"——回测都是上帝视角,前向才是真考验。失败的根因是外生预测(GDP/人口)被 COVID 打废。
8. Ricardian 推广版:把 Pareto Optimal 升级为 Nash Equilibrium(禁止国家间竞争/合谋,点名 2020 俄沙石油战),作者务实承认"收集全部数据要十年",选择"if it's not broken, don't fix it"。

Ricardian 推广版:未尽的野心

Discussion 段末尾,作者给出了模型的"通用版"——Ricardian 模型。当目标国是封闭系统(如 EU 自给自足),naïve 模型够用;但面对美国、加拿大这种农业出口大国,必须扩展成多国贸易系统。

Ricardian 模型来自 David Ricardo 1817 年《On the Principles of Political Economy and Taxation》第 7 章的比较优势理论(英葡葡萄酒-布匹贸易)。在 Smart Farmers 里,不同国家按各自成本和种植面积(比较优势)成比例地生产同种作物,所有国家共享一个国际价格(理想来自 CME/ICE 期货)。

推广版加强了两条假设:

1. 数据可用性假设:Ricardian 模型只用一个国际价格(不分国别),并假设全球统一单产(抹平农技差异)。这要求把运费、关税、配额都算进套利成本——这部分数据"buried inside government documents and updated frequently due to trade frictions",作者吐槽"likely to take a decade for me to collect everything"。
2. Homo Economicus 假设升级:naïve 模型要求 Pareto Optimal,推广版要求 Pareto Optimal 等于 Nash Equilibrium——禁止国家间竞争和合谋。作者点名 2020 俄沙石油战:他们的市场份额竞争导致油价暴跌,合谋又搞垮了 Permian 盆地的页岩油商。这是零和博弈,蛋糕变小。

作者最后很务实:"why bother fixing it when it's not broken?"naïve 模型在 65% 作物上有效,前向预测只要输入(GDP/人口/地/单产)靠谱就能用,没必要为了完美去搞一个十年工程。这是科研的成熟——承认边界,在边界内做到最好。

这种务实也呼应了全章的隐线:Smart Farmers 的价值不在"它是个能赚钱的交易策略",而在"它是个清晰可复现的研究框架"。你拿到这套代码,换成泰国、越南、印尼的数据,换个时间窗,就能复现一遍。这种可复现性,是 quant-trading 全书其他章节都达不到的——它是 Smart Farmers 作为"第二个高潮大项目★"的真正分量。

全章收束

走到这里,我们读完了 Smart Farmers 的全部代码:cleanse data.py 把 FAO 14GB 数据变成 30 类作物的整洁表;country selection.py 用出口占比筛控制组;estimate demand.py 用工具变量 2SLS + constrained QP 估出 α/β/γ;forecast.py 把这些系数塞进 CVXOPT 解出最优产量,再用 Nelder-Mead 找最优成本,最后前向预测 2019-2025 并用棕榈油期货验证。

Smart Farmers 给我们留下三样东西:

1. 一套建模范式:把现实决策(农民种地)翻译成凸优化,价格内生是 QP 的源头,工具变量解决不可观测,合成控制法补缺失值。这套范式可以挪到任何"代表性主体在约束下最优化"的场景。
2. 一个工程模板:CVXOPT 的 P/q/G/h/A/b 标准形式转换、scipy.optimize 与 cvxopt 的嵌套、双轴可视化、外生假设的显式声明。这是 Python 凸优化项目的标准件。
3. 一个思想教训:"I have never seen a bad backtest"。回测是上帝视角,前向才是真考验。下一章 Monte Carlo 批判反思会从另一个角度深化这个教训——所有漂亮的回测都可能只是过拟合。

旅程下一站:第 8 章 Monte Carlo 批判反思与总结。我们会用重采样、白噪声、过拟合诊断等工具,系统性地怀疑前面 7 章所有的"漂亮回测",并给整段旅程画上句号。Smart Farmers 的前向失败会再次出现,作为"回测好 ≠ 策略好"的头号案例。