Notes of EK2002: http://www.humoon.xyz/notes/classic-papers/trade/EK2002/_book/index.html

EK2002 Code (R version)

本项目用 R 语言重写了 EK2002¹ 的祖传 GAUSS 代码，并以文学化编程的方式展示该论文数据处理的全流程。

文学化编程（Literate Programming）：自动生成输出文件（html、pdf、docx 等格式），将文段、图表、代码和代码运行结果混排。或许可以意译为“文档化编程”。

GAUSS 是一门非常古老的数学和统计编程语言，目前已基本无人使用，成为一门僵尸语言。僵尸语言最大的问题还不是技术上的落后，而是帮助文档和资料的极度稀缺，以及生态上的不足。比如，宇宙第一编辑器 VSCode 就不支持 GUASS 语言的代码高亮，每次要看源代码都不得不断网、修改系统时间，才能打开破解版 GAUSS 软件。

EK2002 源代码里的很多命令和函数，我都不知道是什么意思，因为不可能有时间系统地学习一门僵尸语言，或一个一个细查官方文档。更多的是在 GAUSS 软件中调试，逆推那些命令和函数的作用。有时候程序对特定数据的依赖性比较强，无法即时输出，就只能根据运行结果，猜测中间的步骤。以上种种，本项目堪称一个逆向工程。

在宏观层次，逆向工程需要思想的指导，如果我不懂 EK2002 这篇论文建模-估计-模拟的整体架构，和数学上的一系列难点，是绝不可能成功复现其结果的。在微观层次，逆向工程不需要严格一致的还原，而仅需要逻辑等价的还原。EK2002 的源码有很多繁冗、难解之处，如果更简洁的代码能够实现同样的结果，就没有必要 100% 地“翻译”源码。

程序背后涉及的几个理论难点

1. 一般均衡模型的矩阵形式。只有先写成矩阵形式，才能理清逻辑，知道如何进行向量化编程。
2. 根据理论模型设定协方差矩阵的形式，然后根据残差矩阵估计协方差的大小，最后在这个不符合经典 OLS 假设的情况下用 GLS（广义最小二乘法）估计计量模型。
3. 模拟（simulation）的核心是求解非线性方程组，这是数值计算（numerical analysis/computing）领域的典型问题。非线性方程组一般没有解析解，需要用一定的算法寻找数值解。最常用的方法是 Newton-Raphson 迭代法。

源码的不足及解决思路

EK2002 的源码触犯了多条软件工程的禁忌，如：

重写选用 R 的理由

R 是研究场景中最适合处理数据的编程语言，因为

1. 它内置了向量、矩阵和数据框三种数据结构，绝大多数相关运算都可以使用定义得非常简洁、完善的内置函数

2. 它默认支持向量化操作（不像 Python 必须通过第三方库才能引入“广播”功能），而向量化操作可以大大简化向量和矩阵的运算。

   1. 本项目使我深刻认识到矩阵代数的强大。用矩阵代数刻画模型，不仅简洁，而且在逻辑上非常清晰。尤其在国际贸易理论中，涉及多国、多商品、多要素，如果模型写成标量形式，上下标会非常繁琐；写成矩阵形式，变量之间的关系将一目了然。
   2. 对默认支持向量化操作的 R 语言来说，矩阵形式的模型和推导近乎伪代码，稍加改动即称为可以运行的代码。

3. 它对泛函编程范式和管道传输的支持非常好

   1. 数据处理的内在逻辑就是把数据从原始形式经过一系列变换，变为我们想要的形式。在这里，状态性的数据结构、对象都是第二位的，动作性的变换才是第一位的——以至于国外也将数据处理称为 data transformation——所以数据处理天然契合重视过程的泛函编程范式。泛函编程范式的代表，是高阶函数和匿名函数。
   2. 泛函编程范式只有与管道传输联合起来，才能清晰地表现出一系列变换前后接力地作用于数据的链式操作，将数据处理的“流”程表达得清晰而优雅。R 生态中无论第三方模块还是 R base 的最新版本，都对管道传输有良好支持。
   3. R 在面向对象（OOP）的编程范式方面的确比较弱；但在泛函编程范式上，R 不弱于（如果不是更强的话）任何主流编程语言。

4. 无论在计量经济学、还是数值分析（Numerical Analysis，国内也称为数值计算）领域，R 都具有相当完善的生态系统，第三方包提供了许多功能强大、且封装得非常简洁的函数。在执行 EK2002 中的回归和求解非线性方程组等任务时，用 GAUSS 需要写很多行的代码，R 往往三五行即可完成。代码的简洁，能使论文逻辑得到更好的凸显。

5. 综上，在数据科学的研发场景²，R 是目前最高效、最优雅的编程语言。

例：R 数据处理的代码风格

大量使用向量化操作、泛函编程范式和管道传输，尽量避免出现for、while等显式循环。

​x
1
#######################################################
2
## TABLE VI: States of Technology
3
# 由(30)式估计的 S_i 和各种方法估计的 theta，利用(27)式计算 T_i
4
#######################################################
5

6
# %>% 为管道操作符，表示把上一个函数的变换结果作为参数，传递给下一个函数
7

8
table6 <- theta_list %>% # 多个 theta 值，依次传递给高阶函数 map_dfc()
9
  map_dfc(function(theta) {
10
    # 作为 map_dfc() 参数的匿名函数可以对每个 theta 分别计算一个 19 国技术水平列向量
11
    # _dfc 表示将这些向量纵向 (column, c) 堆叠，组成一个 data.frame(df) 类的数据表
12
    # 每行对应一个国家，每列代表一个 theta 情境，这就是 TABLE VI 的雏形
13
    
14
    # 根据 (27) 式 S_i = ln(T_i) / beta - theta * ln(w_i) 求 T_i
15
    tech <- log(adjusted_relative_wage) %>% # ln(w_i)
16
      `*`(theta) %>% # * theta
17
      `+`(source_parameters) %>% # + S_i
18
      `*`(beta) %>% # × beta
19
      exp() # e^
20
    # 上面几行的四则、指数、对数运算全部是向量化操作
21
    # 这样写不仅清楚（每一行包含且仅包含一个运算），而且符合人类面对 (27) 式时的思维：
22
    # 先计算最右边的一坨，移项到等号左边化减为加，两边乘以 beta，最后再做指数变换
23
    
24
    # 如果不用管道，就不得不写成这样
25
    tech <- exp(beta * (source_parameters + theta * log(adjusted_relative_wage)))
26
    # 不仅可读性差，而且只有先在纸上完成对 (27) 的变形，导出 T_i，才能写出这行代码
27
    
28
    
29
    (tech / tech[N]) %>% # 向量化操作，各国技术水平以美国为基准进行标准化
30
      sprintf("%.2f", .) %>%
31
      as.data.frame() %>%
32
      set_colnames(str_c("$\\theta=", sprintf("%.2f", theta), "$"))
33
  }) %>%
34
  mutate( # 再加入两列，一列是国家名称，一列是 S_i (Source-country estimates)
35
    Country = country_table$name,
36
    `Estimated Source-country Competitiveness` = source_parameters %>% sprintf("%.2f", .)
37
  ) %>%
38
  select(Country, `Estimated Source-country Competitiveness`, everything()) 
39
  # 把新加的两列重新排列在最左边

尚未解决的问题

1. 作者在 source code 中用 ${\beta }^{\beta }(1-\beta {)}^{1-\beta }$$\beta^\beta(1-\beta)^{1-\beta}$ 作为 $\gamma$$\gamma$ 值，理由何在？(9) 式中，$\gamma$$\gamma$ 明明是由 $\theta$$\theta$ 和 $\sigma$$\sigma$ 决定的。
2. 劳动力不可跨部门流动的情境中，模拟结果与原文结果的差别很小；而劳动力可以跨部门流动的情境中，模拟结果与原文的差别稍大。何以如此？immobile labor 情境的代码还是基于 mobile labor 情境的，如果后者方向错了，前者按理说也应该是错的；如果前者是对的，后者更应该是对的。
3. Section 6.4 加入关税后的一般均衡模型，正文中没有给出推导。我在 ./docs/ 文件夹下推了一部分，不知正确否。从 GAUSS source code 中应该能完整地逆推出来，较繁琐，待完成。