EK2002: Estimation
source("../config/Rmarkdown_config.R") # 依赖包的导入和一系列 Rmarkdown 设置EK2002 作者使用的编程软件是 GAUSS，其数据文件大都是二进制的，无法用文本编辑器打开，只能用 GAUSS 打印到屏幕上，然后手动复制粘贴到文本文件中，以获得本项目的数据源。
然而，图形化界面的 GAUSS IDE，只能显示小数点后三位，这意味着手动复制粘贴的数据从一开始就是不精确的，会影响所有结果的精度。
只有在命令行风格的 GAUSS Terminal Interface 中打印，才能显示数据的完整位数。从这里复制粘贴，才能保证本项目中数据的精确度。
# 一些标量常数
scalar <- jsonlite::fromJSON("../data/scalar.json")
# 国家代码表
country_table <- fread("../data/country_code.csv")
# 用于参数估计的数据，数值型变量大多经过了对数标准化
regression_data <- fread("../data/regression_data.csv")
# 19国50产业价格数据
price_table <- fread("../data/price.csv")看一下数据结构：
country_table %>% prettify()| id | code | name |
|---|---|---|
| 1 | AL | Australia |
| 2 | AU | Austria |
| 3 | BE | Belgium |
| 4 | CA | Canada |
| 5 | DK | Denmark |
| 6 | FI | Finland |
| 7 | FR | France |
| 8 | GE | Germany |
| 9 | GR | Greece |
| 10 | IT | Italy |
| 11 | JP | Japan |
| 12 | NE | Netherlands |
| 13 | NZ | New Zealand |
| 14 | NO | Norway |
| 15 | PO | Portugal |
| 16 | SP | Spain |
| 17 | SW | Sweden |
| 18 | UK | United Kingdom |
| 19 | US | United States |
regression_dataprice_table## scalar.json
N <- scalar$N # 国家数量
beta <- scalar$beta # 各国平均 beta
theta_estimates <- scalar$theta # 三种不同水平的theta值## regression_data.csv
# 进出口国 index
c(import_country, export_country) %<-% regression_data[, c(1, 2)]
# 对数标准化的制造业双边贸易数据 ln(Xni/Xnn)
trade <- regression_data[, 4] %>% pull()
# ln(X'ni/X'nn), (26)式等号的左边
trade_prime <- regression_data[, 5] %>% pull()
# 由ln(X'ni)的定义可计算(12)式的左边 ln (Xni/Xn)/(Xii/Xi)
# 可将其视为 country n's normalized import share from country i
normalized_trade_share <- -(trade_prime - trade) * beta / (1 - beta) + trade
# 六档距离，说明见P21
c(dist1, dist2, dist3, dist4, dist5, dist6) %<-% regression_data[, 6:11]
# 是否share border/language, 是否均在 EEC/EFTA 中
# 欧共体EEC是EU的前身，EFTA是欧自联
c(border, language, EU, EFTA) %<-% regression_data[, c(12, 18, 19, 20)]
# 横向合并distance, border, language, RTA 等虚拟变量
# 注意，六档距离只需要5个dummy variable
geography_dummy <- cbind(dist2, dist3, dist4, dist5, dist6, border, language, EU, EFTA)
# 以下5个变量的数值都是对数标准化的，即 ln(var_i/var_n)
# 分别是(1)R&D 支出, (2)以平均受教育年限来衡量的 human capital
# (3)人口密度, (4)制造业劳动力数量, (5)工资（均以美元计）
c(r_and_d, human_capital, density, labor, wage) %<-% regression_data[, 13:17]
# 国家首都之间的距离
distance <- regression_data$distance
ln_distance <- log(distance) # 化为对数距离
# 进口国（目的地国）虚拟变量的稀疏矩阵
destination_matrix <- kronecker(diag(N), rep(1, N)) # 克罗内克积
destination_matrix %>%
as.data.table() %>%
set_colnames(1:19 %>% str_c("m", .))# 出口国（来源国）虚拟变量的稀疏矩阵
source_matrix <- kronecker(rep(1, N), diag(N))
source_matrix %>%
as.data.table() %>%
set_colnames(1:19 %>% str_c("S", .))# 来源国虚拟变量 S_i-S_n, (28)式
source_dummy <- source_matrix - destination_matrix
# 这个减法决定了 source_dummy 和 destination_dummy 在(30)式中是不对称的
source_dummy %>% as.data.table() # 观察一下全貌# 目的地国虚拟变量 m_n, (29)式
destination_dummy <- destination_matrix
# 虚拟变量多一个自由度，要以其中一个为基准
# 用其他虚拟变量与这个基准的相对差异，作为回归模型的自变量
# 以美国（第19列）为基准，前18列都减去19列
relative_source_dummy <- (source_dummy[, -N] - source_dummy[, N]) %>%
as.data.table() %>% # 矩阵转换为DT
set_colnames(1:18 %>% str_c("S", .))
relative_source_dummy # 观察一下全貌relative_destination_dummy <- (destination_dummy[, -N] - destination_dummy[, N]) %>%
as.data.table() %>%
set_colnames(1:18 %>% str_c("m", .))
relative_destination_dummy # 观察一下全貌## price.csv
# 50种商品的对数标准化价格（以美国为基准）
price <- price_table %>%
select(-country) %>%
as.matrix()
# 目的地国50种商品对数标准化价格 ln p_n(j)
destination_price <- kronecker(price, rep(1, N))
# 来源地国50种商品对数标准化价格 ln p_i(j)
source_price <- kronecker(rep(1, N), price)
# r_{ni}(j)=ln p_n(j)-ln p_i(j)
relative_price <- destination_price - source_price
# 各行最大的r_{ni}(j)
max_price <- apply(relative_price, 1, max)
# 各行最大r_{ni}(j)的j（商品种类index）
max_price_index <- apply(relative_price, 1, which.max)
# 遍历各行，将最大的r_{ni}(j)替换为-100
for (k in 1:nrow(relative_price)) {
relative_price[k, max_price_index[k]] <- -100
}
# 然后求各行第二大的r_{ni}(j)，作为作为ln d_{ni}的代理变量
ln_dni <- max_price2 <- apply(relative_price, 1, max)
# 价格指数
# 用50种商品对数标准化价格的平均值作为各国价格指数的代理变量
mean_price <- apply(price, 1, mean)
# 目的地国价格指数 ln p_n
destination_price_index <- kronecker(mean_price, rep(1, N))
# 来源地国价格指数 ln p_i
source_price_index <- kronecker(rep(1, N), mean_price)
# (13)式: D_{ni}=ln(p_i)-ln(p_n)+ln(d_{ni})
Dni <- source_price_index - destination_price_index + ln_dni## delete invalid observations
# 删掉 n==i 的行，保留342行
# 因为这些行自变量和因变量的标准化值都等于0，属于无效数据
# 无效行的index
invalid_line_index <- which(import_country == export_country)
# destinations and sources index
import_country_valid <- import_country[-invalid_line_index]
export_country_valid <- export_country[-invalid_line_index]
# ln (Xni/Xn)/(Xii/Xi)
normalized_trade_share_valid <- normalized_trade_share[-invalid_line_index]
# ln(X'ni/X'nn)
trade_prime_valid <- trade_prime[-invalid_line_index]
# D_{ni}
Dni_valid <- Dni[-invalid_line_index]
# exp(D_{ni})，用于 Figure 2
price_measure <- exp(Dni_valid)
# ln_distance
ln_distance_valid <- ln_distance[-invalid_line_index]
# 合并距离、边界、语言、RTA等虚拟变量
geography_dummy_valid <- geography_dummy[-invalid_line_index, ]
# destinations and sources dummies 组合为数据框
relative_source_dummy_valid <- relative_source_dummy[-invalid_line_index, ]
relative_dest_dummy_valid <- relative_destination_dummy[-invalid_line_index, ]
country_dummy_valid <- cbind(
relative_dest_dummy_valid,
relative_source_dummy_valid
) %>% as.data.table()# 标准化贸易份额的最大值
exp(normalized_trade_share_valid) %>% max()#> [1] 0.3584898# 横跨近四个数量级
exp(normalized_trade_share_valid) %>%
max() %>%
`/`(exp(normalized_trade_share_valid) %>% min()) %>%
log() %>%
`/`(log(10))#> [1] 3.985887# 相关系数 correlation
cor(Dni_valid, normalized_trade_share_valid)#> [1] -0.4041783price_measure_table <- # 三列：进口国index，出口国index，exp(Dni)
cbind(import_country_valid, export_country_valid, price_measure) %>%
as.data.table() %>% # 转换数据类型
# 左连接匹配国家代码（code列）
left_join(country_table, by = c("import_country_valid" = "id")) %>%
rename(dest = code) %>%
left_join(country_table, by = c("export_country_valid" = "id")) %>%
rename(source = code) %>%
# 挑选这三列，其余删掉
select(dest, source, price_measure)
price_measure_tabletable2 <- country_table$code %>% # 迭代对象: 19个国家代码
map_dfr( # 每轮迭代生成一行，最后把这些行纵向摞起来
function(country_code) {
# select 进口国为country_code的所有的行
dest_table <- price_measure_table[dest == country_code, ]
source_min <- dest_table[
which.min(price_measure), # select price_measure 最小值所在的行
][
, concatenation := str_c(source, "(", sprintf("%0.2f", price_measure), ")")
] %>% pull(concatenation)
source_max <- dest_table[
which.max(price_measure),
][
, concatenation := str_c(source, "(", sprintf("%0.2f", price_measure), ")")
] %>% pull(concatenation)
# 出口国为country_code的所有的行
source_table <- price_measure_table[source == country_code, ]
dest_min <- source_table[
which.min(price_measure),
][
, concatenation := str_c(dest, "(", sprintf("%0.2f", price_measure), ")")
] %>% pull(concatenation)
dest_max <- source_table[
which.max(price_measure),
][
, concatenation := str_c(dest, "(", sprintf("%0.2f", price_measure), ")")
] %>% pull(concatenation)
# map_dfr()的每个返回值都是Table II的一行
data.table(
Country = str_c(
country_table[code == country_code, ] %>% pull(name),
" (", country_code, ")"
),
source_min, source_max, dest_min, dest_max
)
}
)
table2 %>%
set_colnames(c("Country", "Minimum", "Maximum", "Minimum", "Maximum")) %>%
prettify(caption = "TABLE II \\\n Price Measure Statistics") %>%
add_header_above(c(" " = 1, "Foreign Source \n (Column 1 As Destination)" = 2, "Foreign Destinations \n (Column 1 As Source)" = 2)) %>%
footnote(
general = "本表中的数值为 $e^{D_{ni}}$，其中 $D_{ni}$ 在 (13) 式中定义，为 $\\ln (p_i d_{ni}/p_n)$ 的代理变量。因此，$e^{D_{ni}}$ 表示进口国 $n$ 若全部从出口国 $i$ 的市场购买商品（花费当地价格及运输成本），比实际的购买（在本国市场充分遵循比较优势）会贵多少。",
general_title = "\n Note: ",
footnote_as_chunk = T
)| Country | Minimum | Maximum | Minimum | Maximum |
|---|---|---|---|---|
| Australia (AL) | NE(1.44) | PO(2.25) | BE(1.41) | US(2.03) |
| Austria (AU) | SW(1.39) | NZ(2.16) | UK(1.47) | JP(1.97) |
| Belgium (BE) | GE(1.25) | JP(2.02) | GE(1.35) | SW(1.77) |
| Canada (CA) | US(1.58) | NZ(2.57) | AU(1.57) | US(2.14) |
| Denmark (DK) | FI(1.36) | PO(2.21) | NE(1.48) | US(2.41) |
| Finland (FI) | SW(1.38) | PO(2.61) | DK(1.36) | US(2.87) |
| France (FR) | GE(1.33) | NZ(2.42) | BE(1.40) | JP(2.40) |
| Germany (GE) | BE(1.35) | NZ(2.28) | BE(1.25) | US(2.22) |
| Greece (GR) | SP(1.61) | NZ(2.71) | NE(1.48) | US(2.27) |
| Italy (IT) | FR(1.45) | NZ(2.19) | AU(1.46) | JP(2.10) |
| Japan (JP) | BE(1.62) | PO(3.25) | AL(1.72) | US(3.08) |
| Netherlands (NE) | GE(1.30) | NZ(2.17) | DK(1.39) | NZ(2.01) |
| New Zealand (NZ) | CA(1.60) | PO(2.08) | AL(1.64) | GR(2.71) |
| Norway (NO) | FI(1.45) | JP(2.84) | SW(1.36) | US(2.31) |
| Portugal (PO) | BE(1.49) | JP(2.56) | SP(1.59) | JP(3.25) |
| Spain (SP) | BE(1.39) | JP(2.47) | NO(1.51) | JP(3.05) |
| Sweden (SW) | NO(1.36) | US(2.70) | FI(1.38) | US(2.01) |
| United Kingdom (UK) | NE(1.46) | JP(2.37) | FR(1.52) | NZ(2.04) |
| United States (US) | FR(1.57) | JP(3.08) | CA(1.58) | SW(2.70) |
|
Note: 本表中的数值为 \(e^{D_{ni}}\)，其中 \(D_{ni}\) 在 (13) 式中定义，为 \(\ln (p_i d_{ni}/p_n)\) 的代理变量。因此，\(e^{D_{ni}}\) 表示进口国 \(n\) 若全部从出口国 \(i\) 的市场购买商品（花费当地价格及运输成本），比实际的购买（在本国市场充分遵循比较优势）会贵多少。 |
以 \(\ln (Xni/Xn)/(Xii/Xi)\) 为被解释变量，以 \(D_{ni}\) 为解释变量，估计出来的系数是 \(-\theta\)，因此 \(\theta\) 值还要取一个负号。
最简单的一阶矩估计
x_center <- mean(Dni_valid)
y_center <- mean(normalized_trade_share_valid)
theta <- y_center / x_center
theta#> [1] -8.27593\(\theta=8.28\)
# 无截距项
lm_1 <- lm(normalized_trade_share_valid ~ 0 + Dni_valid)
summary(lm_1) %>%
`[[`("coefficients") %>%
prettify()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| Dni_valid | -8.025323 | 0.1542 | -52.04491 | 0 |
\(\theta=8.03(0.15)\), 与矩估计的估计值很接近。
# 有截距项
lm_2 <- lm(normalized_trade_share_valid ~ Dni_valid)
summary(lm_2) %>%
`[[`("coefficients") %>%
prettify()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -2.173133 | 0.3409387 | -6.373971 | 0 |
| Dni_valid | -4.571893 | 0.5611163 | -8.147851 | 0 |
\(\theta=4.57(0.56)\). 截距项为负，与理论预测（为0）不一致，表明数据中可能有较大的误差.
探讨 \(D_{ni}\) 三个组成部分 \(\ln p_i, \ln p_n, \ln d_{ni}\) 对标准化贸易份额的分别贡献。
# 三个组成部分分别为三列，去掉无效行
trade_share_components <- data.table(
y = normalized_trade_share_valid,
ln_pi = source_price_index[-invalid_line_index],
ln_pn = destination_price_index[-invalid_line_index],
ln_dni = ln_dni[-invalid_line_index]
)
# 无截距项
lm_3 <- lm(y ~ 0 + ln_pi + ln_pn + ln_dni, data = trade_share_components)
summary(lm_3) %>%
`[[`("coefficients") %>%
prettify()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| ln_pi | -9.044924 | 0.4755234 | -19.020987 | 0 |
| ln_pn | 6.350751 | 0.8031414 | 7.907389 | 0 |
| ln_dni | -6.837406 | 0.3272661 | -20.892493 | 0 |
# 有截距项
lm_4 <- lm(y ~ ln_pi + ln_pn + ln_dni, data = trade_share_components)
summary(lm_4) %>%
`[[`("coefficients") %>%
prettify()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -2.358243 | 0.5254622 | -4.487940 | 9.9e-06 |
| ln_pi | -4.878989 | 1.0371534 | -4.704211 | 3.7e-06 |
| ln_pn | 5.488339 | 0.8046659 | 6.820644 | 0.0e+00 |
| ln_dni | -4.555882 | 0.5998471 | -7.595072 | 0.0e+00 |
无论有无截距项，估计出来的参数的正负与理论预测的一致
有截距 OLS 回归的截距不为0，表明我们的价格数据可能有系统误差。
为 price measure \(D_{ni}\) 寻找工具变量（IV），用工具变量法估计 \(\theta\) (以下两种估计仅在源代码中出现，未在论文中出现).
# ivreg::ivreg() 可以做工具变量回归，第一个参数（字符串）为回归公式，形式为：
# 被解释变量 ~ 内生变量 + 外生变量 | 工具变量 + 外生变量
fit_iv_1 <- ivreg(normalized_trade_share_valid ~ Dni_valid | ln_distance_valid)
summary(fit_iv_1, test = TRUE) # 对IV回归进行诊断#>
#> Call:
#> ivreg(formula = normalized_trade_share_valid ~ Dni_valid | ln_distance_valid)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -6.2553 -1.6047 -0.2018 1.1491 9.1736
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3.6579 0.9311 3.929 0.000104 ***
#> Dni_valid -14.5108 1.5732 -9.224 < 2e-16 ***
#>
#> Diagnostic tests:
#> df1 df2 statistic p-value
#> Weak instruments 1 340 110.1 <2e-16 ***
#> Wu-Hausman 1 339 144.5 <2e-16 ***
#> Sargan 0 NA NA NA
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.274 on 340 degrees of freedom
#> Multiple R-Squared: -0.6087, Adjusted R-squared: -0.6134
#> Wald test: 85.08 on 1 and 340 DF, p-value: < 2.2e-16IV 回归通过了诊断，\(\theta=14.51(1.57)\)
instrument_data <- cbind(ln_distance_valid, country_dummy_valid)
formula_instrument <- colnames(instrument_data) %>% str_c(collapse = " + ")
formula_instrument#> [1] "ln_distance_valid + m1 + m2 + m3 + m4 + m5 + m6 + m7 + m8 + m9 + m10 + m11 + m12 + m13 + m14 + m15 + m16 + m17 + m18 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11 + S12 + S13 + S14 + S15 + S16 + S17 + S18"fit_iv_2 <- ivreg(
formula = str_c("normalized_trade_share_valid ~ Dni_valid | ", formula_instrument),
data = cbind(normalized_trade_share_valid, Dni_valid, instrument_data)
)
summary(fit_iv_2, test = TRUE)#>
#> Call:
#> ivreg(formula = str_c("normalized_trade_share_valid ~ Dni_valid | ",
#> formula_instrument), data = cbind(normalized_trade_share_valid,
#> Dni_valid, instrument_data))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -5.19771 -1.18245 -0.03927 1.19787 5.08399
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.8079 0.4281 -1.887 0.06 .
#> Dni_valid -6.8989 0.7130 -9.676 <2e-16 ***
#>
#> Diagnostic tests:
#> df1 df2 statistic p-value
#> Weak instruments 37 304 15.30 < 2e-16 ***
#> Wu-Hausman 1 339 35.26 7.15e-09 ***
#> Sargan 36 NA 263.45 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.681 on 340 degrees of freedom
#> Multiple R-Squared: 0.121, Adjusted R-squared: 0.1185
#> Wald test: 93.62 on 1 and 340 DF, p-value: < 2.2e-16\(\theta=6.90(0.71)\)
估计一个结果被报告在正文中、但方程从未出现过的计量模型2：以 \(\ln (X'_{ni}/X'_{nn})\) 为被解释变量，\(D_{ni}\) 为内生解释变量，source & destination dummies 为外生解释变量
lm_5 <- lm(
trade_prime_valid ~ ., # . 代表所有其他变量
data = cbind(trade_prime_valid, Dni_valid, country_dummy_valid)
)
summary(lm_5) %>%
`[[`("coefficients") %>%
prettify()| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | -3.4199162 | 0.2933648 | -11.6575540 | 0.0000000 |
| Dni_valid | -2.4467276 | 0.4931804 | -4.9611210 | 0.0000012 |
| m1 | -3.0808876 | 0.3050435 | -10.0998319 | 0.0000000 |
| m2 | -0.9392800 | 0.3064023 | -3.0655128 | 0.0023680 |
| m3 | 1.7171125 | 0.3210312 | 5.3487408 | 0.0000002 |
| m4 | -1.4324888 | 0.2989970 | -4.7909809 | 0.0000026 |
| m5 | 0.2877993 | 0.3039180 | 0.9469635 | 0.3444095 |
| m6 | -0.5447744 | 0.2988159 | -1.8231105 | 0.0692687 |
| m7 | 1.2016105 | 0.3008315 | 3.9942971 | 0.0000814 |
| m8 | 2.0422942 | 0.3060856 | 6.6722978 | 0.0000000 |
| m9 | -1.6337823 | 0.3039547 | -5.3750855 | 0.0000002 |
| m10 | 0.4851516 | 0.3055756 | 1.5876645 | 0.1134017 |
| m11 | -0.4603440 | 0.3458418 | -1.3310826 | 0.1841592 |
| m12 | 1.8317076 | 0.3058050 | 5.9897889 | 0.0000000 |
| m13 | -2.2669667 | 0.3286080 | -6.8986952 | 0.0000000 |
| m14 | 0.0388581 | 0.2987892 | 0.1300520 | 0.8966113 |
| m15 | -0.3173925 | 0.3118610 | -1.0177370 | 0.3096122 |
| m16 | -0.7589797 | 0.3012969 | -2.5190428 | 0.0122802 |
| m17 | 1.0137428 | 0.2987878 | 3.3928520 | 0.0007833 |
| m18 | 1.8161109 | 0.3006685 | 6.0402427 | 0.0000000 |
| S1 | -1.4140555 | 0.2066347 | -6.8432638 | 0.0000000 |
| S2 | -0.7551922 | 0.2073354 | -3.6423696 | 0.0003175 |
| S3 | -3.0754706 | 0.2161052 | -14.2313592 | 0.0000000 |
| S4 | -0.7201263 | 0.2058235 | -3.4987570 | 0.0005373 |
| S5 | -1.3617033 | 0.2077610 | -6.5541802 | 0.0000000 |
| S6 | -0.1235467 | 0.2056388 | -0.6007946 | 0.5484244 |
| S7 | 1.7854210 | 0.2061861 | 8.6592711 | 0.0000000 |
| S8 | 2.8383352 | 0.2095910 | 13.5422546 | 0.0000000 |
| S9 | -2.4888271 | 0.2066605 | -12.0430736 | 0.0000000 |
| S10 | 1.9614312 | 0.2090452 | 9.3828110 | 0.0000000 |
| S11 | 3.1925417 | 0.2246246 | 14.2127856 | 0.0000000 |
| S12 | -1.7780664 | 0.2082723 | -8.5372199 | 0.0000000 |
| S13 | -2.1422594 | 0.2347276 | -9.1265782 | 0.0000000 |
| S14 | -0.9624557 | 0.2075509 | -4.6372043 | 0.0000053 |
| S15 | -0.9158288 | 0.2234968 | -4.0977259 | 0.0000536 |
| S16 | 0.5262126 | 0.2061498 | 2.5525743 | 0.0111817 |
| S17 | 0.4185237 | 0.2070864 | 2.0210099 | 0.0441550 |
| S18 | 1.8218525 | 0.2074390 | 8.7825923 | 0.0000000 |
\(\theta=2.45(0.49)\)
以 geography_dummies 为工具变量，包括距离、边界、语言、RTA 等，解释内生变量 D_ni；country_dummies 为外生变量。
TSLS_data <- cbind(trade_prime_valid, Dni_valid, geography_dummy_valid, country_dummy_valid)
TSLS_data# 外生变量
formula_exogenous <- colnames(country_dummy_valid) %>% str_c(collapse = " + ")
formula_exogenous#> [1] "m1 + m2 + m3 + m4 + m5 + m6 + m7 + m8 + m9 + m10 + m11 + m12 + m13 + m14 + m15 + m16 + m17 + m18 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11 + S12 + S13 + S14 + S15 + S16 + S17 + S18"# 工具变量
formula_instrument <- colnames(geography_dummy_valid) %>% str_c(collapse = " + ")
formula_instrument#> [1] "dist2 + dist3 + dist4 + dist5 + dist6 + border + language + EU + EFTA"# 回归公式
formula <- str_c(
"trade_prime_valid ~ Dni_valid + ", formula_exogenous, " | ",
formula_instrument, " + ", formula_exogenous
)
formula#> [1] "trade_prime_valid ~ Dni_valid + m1 + m2 + m3 + m4 + m5 + m6 + m7 + m8 + m9 + m10 + m11 + m12 + m13 + m14 + m15 + m16 + m17 + m18 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11 + S12 + S13 + S14 + S15 + S16 + S17 + S18 | dist2 + dist3 + dist4 + dist5 + dist6 + border + language + EU + EFTA + m1 + m2 + m3 + m4 + m5 + m6 + m7 + m8 + m9 + m10 + m11 + m12 + m13 + m14 + m15 + m16 + m17 + m18 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11 + S12 + S13 + S14 + S15 + S16 + S17 + S18"fit_iv_3 <- ivreg(formula = formula, data = TSLS_data)
summary(fit_iv_3, test = TRUE)#>
#> Call:
#> ivreg(formula = formula, data = TSLS_data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -4.03108 -0.87353 -0.07295 0.91475 5.54823
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.690744 1.021376 2.634 0.008860 **
#> Dni_valid -12.862154 1.736076 -7.409 1.27e-12 ***
#> m1 -4.379061 0.516788 -8.474 1.04e-15 ***
#> m2 -2.373090 0.526655 -4.506 9.44e-06 ***
#> m3 -0.762720 0.625382 -1.220 0.223560
#> m4 -1.195235 0.470969 -2.538 0.011654 *
#> m5 -0.886755 0.508505 -1.744 0.082197 .
#> m6 -0.634263 0.469543 -1.351 0.177761
#> m7 0.461957 0.485230 0.952 0.341834
#> m8 0.638974 0.524368 1.219 0.223956
#> m9 -0.455004 0.508776 -0.894 0.371863
#> m10 -0.867693 0.520670 -1.666 0.096645 .
#> m11 3.217836 0.772093 4.168 4.02e-05 ***
#> m12 0.455939 0.522336 0.873 0.383415
#> m13 0.621791 0.672382 0.925 0.355825
#> m14 0.008904 0.469333 0.019 0.984876
#> m15 1.569687 0.564966 2.778 0.005803 **
#> m16 -1.578771 0.488795 -3.230 0.001374 **
#> m17 0.990810 0.469322 2.111 0.035575 *
#> m18 1.106633 0.483976 2.287 0.022910 *
#> S1 -1.842036 0.330782 -5.569 5.65e-08 ***
#> S2 -1.314255 0.336172 -3.909 0.000114 ***
#> S3 -4.478607 0.398783 -11.231 < 2e-16 ***
#> S4 -0.904355 0.324455 -2.787 0.005649 **
#> S5 -1.987282 0.339413 -5.855 1.24e-08 ***
#> S6 -0.116991 0.323000 -0.362 0.717454
#> S7 1.468290 0.327295 4.486 1.03e-05 ***
#> S8 1.982772 0.353079 5.616 4.42e-08 ***
#> S9 -2.055326 0.330982 -6.210 1.74e-09 ***
#> S10 1.167665 0.349046 3.345 0.000925 ***
#> S11 5.101323 0.453382 11.252 < 2e-16 ***
#> S12 -2.475357 0.343273 -7.211 4.43e-12 ***
#> S13 0.247988 0.512895 0.484 0.629084
#> S14 -1.556102 0.337816 -4.606 6.03e-06 ***
#> S15 0.932966 0.446422 2.090 0.037460 *
#> S16 0.219792 0.327012 0.672 0.502017
#> S17 -0.097743 0.334265 -0.292 0.770172
#> S18 1.245900 0.336963 3.697 0.000258 ***
#>
#> Diagnostic tests:
#> df1 df2 statistic p-value
#> Weak instruments 9 296 8.176 7.44e-11 ***
#> Wu-Hausman 1 303 173.953 < 2e-16 ***
#> Sargan 8 NA 47.419 1.28e-07 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.406 on 304 degrees of freedom
#> Multiple R-Squared: 0.7955, Adjusted R-squared: 0.7706
#> Wald test: 38.05 on 37 and 304 DF, p-value: < 2.2e-16\(\theta=12.86(1.74)\)
估计 (30) 式
# 组织数据框
table3_data <- cbind(
trade_prime, dist1, geography_dummy,
relative_source_dummy, relative_destination_dummy
) %>%
as.data.table() %>%
`[`(-invalid_line_index, )
# OLS
lm_table3 <- lm(trade_prime ~ 0 + ., data = table3_data)
summary(lm_table3)#>
#> Call:
#> lm(formula = trade_prime ~ 0 + ., data = table3_data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.37812 -0.28260 0.00324 0.27485 1.35971
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> dist1 -3.102440 0.154903 -20.028 < 2e-16 ***
#> dist2 -3.665941 0.106969 -34.271 < 2e-16 ***
#> dist3 -4.033395 0.095713 -42.141 < 2e-16 ***
#> dist4 -4.218077 0.153080 -27.555 < 2e-16 ***
#> dist5 -6.064372 0.087005 -69.701 < 2e-16 ***
#> dist6 -6.558933 0.097792 -67.070 < 2e-16 ***
#> border 0.303564 0.137159 2.213 0.027645 *
#> language 0.510080 0.145505 3.506 0.000526 ***
#> EU 0.035912 0.121613 0.295 0.767974
#> EFTA 0.536071 0.182563 2.936 0.003581 **
#> S1 0.192558 0.149741 1.286 0.199471
#> S2 -1.161839 0.128054 -9.073 < 2e-16 ***
#> S3 -3.335625 0.118827 -28.071 < 2e-16 ***
#> S4 0.411561 0.141915 2.900 0.004010 **
#> S5 -1.749550 0.121735 -14.372 < 2e-16 ***
#> S6 -0.522466 0.129777 -4.026 7.22e-05 ***
#> S7 1.281384 0.118952 10.772 < 2e-16 ***
#> S8 2.352748 0.121609 19.347 < 2e-16 ***
#> S9 -2.813521 0.122484 -22.971 < 2e-16 ***
#> S10 1.781510 0.120437 14.792 < 2e-16 ***
#> S11 4.198787 0.132413 31.710 < 2e-16 ***
#> S12 -2.189319 0.120839 -18.118 < 2e-16 ***
#> S13 -1.197687 0.149741 -7.998 2.87e-14 ***
#> S14 -1.346152 0.129424 -10.401 < 2e-16 ***
#> S15 -1.573180 0.126130 -12.473 < 2e-16 ***
#> S16 0.303089 0.121268 2.499 0.012984 *
#> S17 0.009816 0.129869 0.076 0.939800
#> S18 1.374167 0.121654 11.296 < 2e-16 ***
#> m1 0.236221 0.256161 0.922 0.357197
#> m2 -1.678415 0.203764 -8.237 5.75e-15 ***
#> m3 1.120118 0.180023 6.222 1.67e-09 ***
#> m4 0.688596 0.237671 2.897 0.004045 **
#> m5 -0.505889 0.187633 -2.696 0.007416 **
#> m6 -1.334488 0.207808 -6.422 5.36e-10 ***
#> m7 0.218294 0.180352 1.210 0.227099
#> m8 0.998812 0.187305 5.333 1.93e-07 ***
#> m9 -2.356411 0.189573 -12.430 < 2e-16 ***
#> m10 0.070178 0.184252 0.381 0.703566
#> m11 1.584890 0.214620 7.385 1.56e-12 ***
#> m12 1.004782 0.185301 5.422 1.22e-07 ***
#> m13 0.066574 0.256161 0.260 0.795129
#> m14 -1.000409 0.207199 -4.828 2.21e-06 ***
#> m15 -1.206780 0.198907 -6.067 3.97e-09 ***
#> m16 -1.156611 0.186420 -6.204 1.85e-09 ***
#> m17 -0.024864 0.208183 -0.119 0.905012
#> m18 0.816808 0.187422 4.358 1.81e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.4913 on 296 degrees of freedom
#> Multiple R-squared: 0.9935, Adjusted R-squared: 0.9925
#> F-statistic: 984.2 on 46 and 296 DF, p-value: < 2.2e-16由理论推断，(30) 式的回归模型违反了 OLS 的基本假设: \(\text{cov}(\varepsilon_j, \varepsilon_k)=E(\varepsilon_{j}\varepsilon_{k})-E(\varepsilon_{j})E(\varepsilon_{k})=E(\varepsilon_{j}\varepsilon_{k}) = 0, \forall j \neq k\).
这是因为，贸易具有互惠性，如果 \(i\) 国向 \(n\) 国出口的贸易障碍偏低（实际值低于模型拟合值），则一般 \(n\) 国向 \(i\) 国出口的贸易障碍也会偏低。这两个误差 \(\varepsilon_{n \leftarrow i}, \varepsilon_{i\leftarrow n}\)（以下简写为 \(\varepsilon_{ni}\) 和 \(\varepsilon_{in}\)）往往同时偏高或偏低，有一定相关性，故 \(\text{cov}(\varepsilon_{ni}, \varepsilon_{in})>0\)
设误差 \(\varepsilon_{ni}=\varepsilon_{ni}^1+\varepsilon_{ni}^2\) 可以分为两部分，两部分之间不相关
第一部分 \(\varepsilon_{ni}^1\) 符合 OLS 的经典假设，即 \(\begin{aligned}E(\varepsilon_{ab}^1\varepsilon_{cd}^1) = \left\{\begin{matrix} \sigma_1^2, & \forall a=c \land b=d\\ 0, & \text{otherwise.} \end{matrix}\right.\end{aligned}\)
第二部分专门体现两国贸易的互惠性，有性质 \(\varepsilon_{ni}^2=\varepsilon_{in}^2\)，使得 \(E(\varepsilon_{ni}^2\varepsilon_{in}^2)= E(\varepsilon_{ni}^2\varepsilon_{ni}^2) =\sigma_2^2\)
则总误差满足 \[\begin{align} E(\varepsilon_{ab}\varepsilon_{cd}) & = E(\varepsilon_{ab}^1\varepsilon_{cd}^1)+E(\varepsilon_{ab}^1\varepsilon_{cd}^2)+E(\varepsilon_{ab}^2\varepsilon_{cd}^1)+E(\varepsilon_{ab}^2\varepsilon_{cd}^2) \\ &=E(\varepsilon_{ab}^1\varepsilon_{cd}^1)+E(\varepsilon_{ab}^2\varepsilon_{cd}^2) \\ &= \left\{\begin{matrix} \sigma_1^2+\sigma_2^2, & \forall a=c \land b=d \\ \sigma_2^2, & \forall a=d \land b=c\\ 0, & \text{otherwise.} \end{matrix}\right. \end{align}\]
误差项有了上述性质，便不能再使用 OLS estimation，而应使用广义最小二乘法（GLS）3。
对于纯粹的异方差问题（误差的协方差矩阵只是对角元素不同，但非对角元素保证为
0），R base 中有 lm() 或 lm.gls()
函数可以直接用；但对于协方差矩阵非对角元素不为零的情况，目前似乎没有很好的
library 提供现成的函数，需要自己写矩阵运算。
# 由 OLS 求得所有观测的残差，并将其扩展为矩阵，作为估计误差协方差的材料
ee_matrix <- lm_table3$residuals %o% lm_table3$residuals
# 这是一个 342*342 的大矩阵
# 估计 sigma_1^2 + sigma_2^2
sum_sigma_square <- diag(ee_matrix) %>% # 从n到i与从n到i的trade_prime的残差之积
mean() # 期望为 sigma_1^2 + sigma_2^2
# 估计 sigma_2^2
observation_valid <- data.table(n = import_country, i = export_country) %>%
`[`(-invalid_line_index, ) %>%
`[`(, row_index := .I)
print(observation_valid)#> n i row_index
#> 1: 1 2 1
#> 2: 1 3 2
#> 3: 1 4 3
#> 4: 1 5 4
#> 5: 1 6 5
#> ---
#> 338: 19 14 338
#> 339: 19 15 339
#> 340: 19 16 340
#> 341: 19 17 341
#> 342: 19 18 342sigma2_square <- observation_valid$row_index %>%
map_dbl(function(row_index) {
row_n <- observation_valid$n[row_index]
row_i <- observation_valid$i[row_index]
col_index <- observation_valid[n == row_i & i == row_n, ]$row_index
ee_matrix[row_index, col_index] # 从n到i与从i到n的trade_prime的残差之积
}) %>%
mean() # 期望为 sigma_2^2
sigma2_square#> [1] 0.05064993sum_sigma_square - sigma2_square#> [1] 0.1582666因此
\(\theta^2\sigma_2^2=\) 0.05
\(\theta^2\sigma_1^2=\) 0.16
# 误差协方差矩阵 Omega
Omega <- diag(nrow(observation_valid)) * sum_sigma_square
observation_valid$row_index %>%
walk(function(row_index) {
row_n <- observation_valid$n[row_index]
row_i <- observation_valid$i[row_index]
col_index <- observation_valid[n == row_i & i == row_n, ]$row_index
Omega[row_index, col_index] <<- sigma2_square # "<<-" 对环境外变量赋值
})
Omega %>% round(2) %>% as.data.table()# 求逆
Omega_inverse <- solve(Omega)
# GLS estimation
Y <- table3_data$trade_prime
X <- table3_data[, -1] %>% as.matrix()
estimate <- solve(t(X) %*% Omega_inverse %*% X) %*% t(X) %*% Omega_inverse %*% Y
estimate_S19 <- -estimate[11:28] %>% sum() # 保证虚拟变量效应之和为0
estimate_m19 <- -estimate[29:46] %>% sum()
var_estimate <- solve(t(X) %*% Omega_inverse %*% X)
stand_error <- diag(var_estimate) %>% sqrt()
stand_error_S19 <- var_estimate[11:28, 11:28] %>%
sum() %>% # -S19等于S1到S18的和，因此S19的方差应等于S1到S18所有的方差和协方差之和
sqrt()
stand_error_m19 <- var_estimate[29:46, 29:46] %>%
sum() %>%
sqrt()# Table III
gls_estimate <- c(
estimate[1:28], estimate_S19,
estimate[29:46], estimate_m19
)
gls_stand_error <- c(
stand_error[1:28], stand_error_S19,
stand_error[29:46], stand_error_m19
)
gls_parameters <- c(
1:6 %>% str_c("$-\\theta d_{", ., "}$"),
"$-\\theta b$", "$-\\theta l$", "$-\\theta e_1$", "$-\\theta e_2$",
1:19 %>% str_c("$S_{", ., "}$"),
1:19 %>% str_c("$-\\theta m_{", ., "}$")
)
variable_names <- c(
"Distance [0, 375)", "Distance [375, 750)", "Distance [750, 1500)",
"Distance [1500, 3000)", "Distance [3000, 6000)", "Distance [6000, maximum)",
"Shared border", "Shared language", "Both in EEC/EU", "Both in EFTA",
country_table$name %>% str_c("Source country: ", .),
country_table$name %>% str_c("Destination country: ", .)
)
table3 <- data.table(
variables = variable_names,
parameters = gls_parameters,
estimate = round(gls_estimate, 2),
std_error = round(gls_stand_error, 2)
)
table3 %>%
mutate(std_error = str_c("(", std_error, ")")) %>%
set_colnames(c("Variable", "parameter", "est.", "s.e.")) %>%
prettify(caption = "TABLE III \\\n Bilateral Trade Equation") %>%
footnote(
general = "Estimated by GLS. The parameter are normalized so that ${\\textstyle \\sum_{i=1}^{19} S_i}=0$ and ${\\textstyle \\sum_{n=1}^{19} m_n}=0$. Standard errors are in parentheses.",
general_title = "\n Note: ",
footnote_as_chunk = T
)| Variable | parameter | est. | s.e. |
|---|---|---|---|
| Distance [0, 375) | \(-\theta d_{1}\) | -3.10 | (0.16) |
| Distance [375, 750) | \(-\theta d_{2}\) | -3.66 | (0.11) |
| Distance [750, 1500) | \(-\theta d_{3}\) | -4.03 | (0.1) |
| Distance [1500, 3000) | \(-\theta d_{4}\) | -4.22 | (0.16) |
| Distance [3000, 6000) | \(-\theta d_{5}\) | -6.06 | (0.09) |
| Distance [6000, maximum) | \(-\theta d_{6}\) | -6.56 | (0.1) |
| Shared border | \(-\theta b\) | 0.30 | (0.14) |
| Shared language | \(-\theta l\) | 0.51 | (0.15) |
| Both in EEC/EU | \(-\theta e_1\) | 0.04 | (0.13) |
| Both in EFTA | \(-\theta e_2\) | 0.54 | (0.19) |
| Source country: Australia | \(S_{1}\) | 0.19 | (0.15) |
| Source country: Austria | \(S_{2}\) | -1.16 | (0.12) |
| Source country: Belgium | \(S_{3}\) | -3.34 | (0.11) |
| Source country: Canada | \(S_{4}\) | 0.41 | (0.14) |
| Source country: Denmark | \(S_{5}\) | -1.75 | (0.12) |
| Source country: Finland | \(S_{6}\) | -0.52 | (0.12) |
| Source country: France | \(S_{7}\) | 1.28 | (0.11) |
| Source country: Germany | \(S_{8}\) | 2.35 | (0.12) |
| Source country: Greece | \(S_{9}\) | -2.81 | (0.12) |
| Source country: Italy | \(S_{10}\) | 1.78 | (0.11) |
| Source country: Japan | \(S_{11}\) | 4.20 | (0.13) |
| Source country: Netherlands | \(S_{12}\) | -2.19 | (0.11) |
| Source country: New Zealand | \(S_{13}\) | -1.20 | (0.15) |
| Source country: Norway | \(S_{14}\) | -1.35 | (0.12) |
| Source country: Portugal | \(S_{15}\) | -1.57 | (0.12) |
| Source country: Spain | \(S_{16}\) | 0.30 | (0.12) |
| Source country: Sweden | \(S_{17}\) | 0.01 | (0.12) |
| Source country: United Kingdom | \(S_{18}\) | 1.37 | (0.12) |
| Source country: United States | \(S_{19}\) | 3.98 | (0.14) |
| Destination country: Australia | \(-\theta m_{1}\) | 0.24 | (0.27) |
| Destination country: Austria | \(-\theta m_{2}\) | -1.68 | (0.21) |
| Destination country: Belgium | \(-\theta m_{3}\) | 1.12 | (0.19) |
| Destination country: Canada | \(-\theta m_{4}\) | 0.69 | (0.25) |
| Destination country: Denmark | \(-\theta m_{5}\) | -0.51 | (0.19) |
| Destination country: Finland | \(-\theta m_{6}\) | -1.33 | (0.22) |
| Destination country: France | \(-\theta m_{7}\) | 0.22 | (0.19) |
| Destination country: Germany | \(-\theta m_{8}\) | 1.00 | (0.19) |
| Destination country: Greece | \(-\theta m_{9}\) | -2.36 | (0.2) |
| Destination country: Italy | \(-\theta m_{10}\) | 0.07 | (0.19) |
| Destination country: Japan | \(-\theta m_{11}\) | 1.59 | (0.22) |
| Destination country: Netherlands | \(-\theta m_{12}\) | 1.00 | (0.19) |
| Destination country: New Zealand | \(-\theta m_{13}\) | 0.07 | (0.27) |
| Destination country: Norway | \(-\theta m_{14}\) | -1.00 | (0.21) |
| Destination country: Portugal | \(-\theta m_{15}\) | -1.21 | (0.21) |
| Destination country: Spain | \(-\theta m_{16}\) | -1.16 | (0.19) |
| Destination country: Sweden | \(-\theta m_{17}\) | -0.02 | (0.22) |
| Destination country: United Kingdom | \(-\theta m_{18}\) | 0.81 | (0.19) |
| Destination country: United States | \(-\theta m_{19}\) | 2.46 | (0.25) |
|
Note: Estimated by GLS. The parameter are normalized so that \({\textstyle \sum_{i=1}^{19} S_i}=0\) and \({\textstyle \sum_{n=1}^{19} m_n}=0\). Standard errors are in parentheses. |
TABLE III 的估计值对于下一步的 simulation 非常重要。
table3_estimate <- gls_estimate # 是用高精度估计值
# table3$estimate 的数据只有两位小数，不能用
barrier_dummy <- cbind(regression_data[, c(6:12, 18:20)], destination_dummy) %>%
as.matrix() # 361*29
barrier_estimate <- gls_estimate[c(1:10, 30:48)] # 29*1
barrier_measure <- barrier_dummy %*% barrier_estimate %>% # 361*1
`[<-`(invalid_line_index, 0) %>% # n==i 的行取0
matrix(nrow = N, byrow = TRUE)# 保存为 R 二进制文件 .rda，使用时直接 load()
save(table3_estimate, barrier_measure,
file = "../data/estimation-output.rda"
)# 所有国家相对于美国的 R&D stock，TABLE 4 的第一列
relative_r_and_d <- exp(r_and_d[1:N]) / exp(r_and_d[N])
# 各国平均教育年限，TABLE 4 的第二列
edu_year <- regression_data[1:N, 22][[1]] # 一定要转换为向量
# 各国相对于美国的有效制造业劳动力，TABLE 4 的第三列
adjusted_relative_labor <- exp(labor[1:N]) * exp(0.06 * edu_year) /
(exp(labor[N]) * exp(0.06 * edu_year[N]))
# 各国相对于美国的有效（经人力资本调整）工资
adjusted_relative_wage <- exp(wage[1:N]) * exp(-0.06 * edu_year) /
(exp(wage[N]) * exp(-0.06 * edu_year[N]))
# 各国相对于美国的人口密度
relative_density <- exp(density[1:N]) / exp(density[N])
table4 <- data.table(
country = country_table$name,
research_stock = relative_r_and_d %>% sprintf("%.4f", .),
years_of_schooling = edu_year %>% sprintf("%.2f", .),
labor_force_HK_adjusted = adjusted_relative_labor %>% sprintf("%.3f", .),
density = relative_density %>% sprintf("%.2f", .)
)
table4 %>%
set_colnames(c(
"Country", "Research \\\n Stock \\\n (U.S.=1)", "Years of\\\n Schooling \\\n (years/person)", "Labor Force \\\n (HK adjusted) \\\n (U.S.=1)", "Density \\\n (pop/area) \\\n (U.S.=1)"
)) %>%
prettify(caption = "TABLE IV \\\n Data for Alternative Parameters") %>%
footnote(
general = "Research stocks, in 1990, are from Coe and Helpman (1995). Average years of schooling $H_i$, in 1985, are from Kyriacou (1991). Labor forces, in 1990, are from Summers and Heston (1991). They are adjusted for human capital by multiplying the country $i$ figure by $e^{0.06H_i}$. See the Appendix for complete definitions.",
general_title = "\n Note: ",
footnote_as_chunk = T
)| Country |
Research Stock (U.S.=1) |
Years of Schooling (years/person) |
Labor Force (HK adjusted) (U.S.=1) |
Density (pop/area) (U.S.=1) |
|---|---|---|---|---|
| Australia | 0.0087 | 8.72 | 0.054 | 0.08 |
| Austria | 0.0063 | 8.57 | 0.024 | 3.43 |
| Belgium | 0.0151 | 9.35 | 0.029 | 12.03 |
| Canada | 0.0300 | 9.98 | 0.095 | 0.10 |
| Denmark | 0.0051 | 6.91 | 0.017 | 4.47 |
| Finland | 0.0053 | 10.83 | 0.019 | 0.55 |
| France | 0.1108 | 9.54 | 0.181 | 3.88 |
| Germany | 0.1683 | 10.33 | 0.225 | 9.50 |
| Greece | 0.0005 | 8.40 | 0.025 | 2.87 |
| Italy | 0.0446 | 9.13 | 0.159 | 7.16 |
| Japan | 0.2493 | 9.47 | 0.544 | 12.42 |
| Netherlands | 0.0278 | 9.47 | 0.043 | 13.64 |
| New Zealand | 0.0010 | 9.27 | 0.010 | 0.47 |
| Norway | 0.0057 | 9.24 | 0.015 | 0.49 |
| Portugal | 0.0007 | 6.52 | 0.026 | 4.01 |
| Spain | 0.0084 | 9.70 | 0.100 | 2.88 |
| Sweden | 0.0207 | 9.64 | 0.031 | 0.71 |
| United Kingdom | 0.1423 | 8.50 | 0.186 | 8.76 |
| United States | 1.0000 | 12.09 | 1.000 | 1.00 |
|
Note: Research stocks, in 1990, are from Coe and Helpman (1995). Average years of schooling \(H_i\), in 1985, are from Kyriacou (1991). Labor forces, in 1990, are from Summers and Heston (1991). They are adjusted for human capital by multiplying the country \(i\) figure by \(e^{0.06H_i}\). See the Appendix for complete definitions. |
估计 Section 5.2 的未编号公式
source_parameters <- gls_estimate[11:29]competitiveness <- data.table(
research_stock = log(relative_r_and_d),
human_capital = 1 / edu_year,
wage = log(adjusted_relative_wage)
) %>%
cbind(source_parameters, .)
lm_table5 <- lm(source_parameters ~ ., data = competitiveness)
summary(lm_table5)#>
#> Call:
#> lm(formula = source_parameters ~ ., data = competitiveness)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.63550 -0.31749 0.09465 0.59003 1.46797
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 5.8130 2.0655 2.814 0.0131 *
#> research_stock 1.0384 0.1728 6.010 2.39e-05 ***
#> human_capital -18.0562 20.5483 -0.879 0.3934
#> wage -2.8472 1.0227 -2.784 0.0139 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.109 on 15 degrees of freedom
#> Multiple R-squared: 0.7703, Adjusted R-squared: 0.7244
#> F-statistic: 16.77 on 3 and 15 DF, p-value: 4.673e-05工资与技术水平之间可能存在内生性问题
将劳动力和人口密度作为工资的工具变量，这二者外生于技术
relative_density <- exp(density[1:N]) / exp(density[N])
competitiveness_2SLS <- data.table(
research_stock = log(relative_r_and_d),
human_capital = 1 / edu_year,
wage = log(adjusted_relative_wage),
labor = log(adjusted_relative_labor),
density = log(relative_density)
) %>%
cbind(source_parameters, .)
fit_iv <- ivreg(
formula = source_parameters ~ research_stock + human_capital + wage |
research_stock + human_capital + labor + density,
data = competitiveness_2SLS
)
summary(fit_iv, test = TRUE)#>
#> Call:
#> ivreg(formula = source_parameters ~ research_stock + human_capital +
#> wage | research_stock + human_capital + labor + density,
#> data = competitiveness_2SLS)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.4724 -0.5927 0.0668 0.6626 1.3383
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.4385 2.1624 2.978 0.00939 **
#> research_stock 1.0941 0.1815 6.027 2.32e-05 ***
#> human_capital -22.6873 21.2475 -1.068 0.30251
#> wage -3.5988 1.2054 -2.986 0.00924 **
#>
#> Diagnostic tests:
#> df1 df2 statistic p-value
#> Weak instruments 2 14 20.538 6.86e-05 ***
#> Wu-Hausman 1 14 1.654 0.2193
#> Sargan 1 NA 3.229 0.0723 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.129 on 15 degrees of freedom
#> Multiple R-Squared: 0.7621, Adjusted R-squared: 0.7145
#> Wald test: 16.67 on 3 and 15 DF, p-value: 4.843e-05估计得 \(\theta=3.60(1.21)\)
由 (30) 式估计的 \(S\) 和各种方法估计的 \(\theta\)，利用 (27) 式计算技术水平 \(T\)（另一种思路是由 Section 5.2 估计的 \(a_0\)、\(a_R\)、\(a_H\) 计算 \(T\)）
table6 <- theta_estimates %>%
map_dfc(function(theta) {
# (27) 式变形
tech <- log(adjusted_relative_wage) %>% # ln(w_i)
`*`(theta) %>% # * theta
`+`(source_parameters) %>% # + S_i
`*`(beta) %>% # × beta
exp() # e^
(tech / tech[N]) %>%
sprintf("%.2f", .) %>%
as.data.frame() %>%
set_colnames(str_c("$\\theta=", sprintf("%.2f", theta), "$"))
}) %>%
mutate(
Country = country_table$name,
`Estimated Source-country Competitiveness` = source_parameters %>% sprintf("%.2f", .)
) %>%
select(Country, `Estimated Source-country Competitiveness`, everything())
table6 %>%
prettify(caption = "TABLE VI \\\n States of Technology") %>%
footnote(
general = "The estimates of source-country competitiveness are the same as those shown in Table III. For an estimated parameter $\\hat{S_i}$, the implied state of technology is $T_i=(e^{\\hat{S_i}}w_i^\\theta)^\\beta$. States of technology are normalized relative to the U.S. value.",
general_title = "\n Note: ",
footnote_as_chunk = T
) %>%
add_header_above(c(
" " = 1, " " = 1,
"Implied States of Technology" = 3
))| Country | Estimated Source-country Competitiveness | \(\theta=3.60\) | \(\theta=8.28\) | \(\theta=12.86\) |
|---|---|---|---|---|
| Australia | 0.19 | 0.36 | 0.27 | 0.20 |
| Austria | -1.16 | 0.30 | 0.26 | 0.23 |
| Belgium | -3.34 | 0.22 | 0.24 | 0.26 |
| Canada | 0.41 | 0.47 | 0.46 | 0.46 |
| Denmark | -1.75 | 0.32 | 0.35 | 0.38 |
| Finland | -0.52 | 0.41 | 0.45 | 0.50 |
| France | 1.28 | 0.60 | 0.64 | 0.69 |
| Germany | 2.35 | 0.75 | 0.81 | 0.86 |
| Greece | -2.81 | 0.14 | 0.07 | 0.04 |
| Italy | 1.78 | 0.57 | 0.50 | 0.45 |
| Japan | 4.20 | 0.97 | 0.89 | 0.81 |
| Netherlands | -2.19 | 0.28 | 0.30 | 0.32 |
| New Zealand | -1.20 | 0.22 | 0.12 | 0.07 |
| Norway | -1.35 | 0.37 | 0.43 | 0.50 |
| Portugal | -1.57 | 0.13 | 0.04 | 0.01 |
| Spain | 0.30 | 0.33 | 0.21 | 0.14 |
| Sweden | 0.01 | 0.47 | 0.51 | 0.57 |
| United Kingdom | 1.37 | 0.53 | 0.49 | 0.44 |
| United States | 3.98 | 1.00 | 1.00 | 1.00 |
|
Note: The estimates of source-country competitiveness are the same as those shown in Table III. For an estimated parameter \(\hat{S_i}\), the implied state of technology is \(T_i=(e^{\hat{S_i}}w_i^\theta)^\beta\). States of technology are normalized relative to the U.S. value. |
用 (29) 式计算贸易障碍
barrier_parameters <- gls_estimate[-11:-29]
barrier_names <- c(
"Distance [0, 375)", "Distance [375, 750)", "Distance [750, 1500)",
"Distance [1500, 3000)", "Distance [3000, 6000)", "Distance [6000, maximum)",
"Shared border", "Shared language", "Both in EEC/EU", "Both in EFTA",
country_table$name
)
table7 <- theta_estimates %>%
map_dfc(function(theta) {
# (29), (30) 变形
(barrier_parameters / (-theta)) %>% # 除以 -theta, 才是对 ln(d_ni) 的影响
exp() %>% # 对dni的影响
`-`(1) %>%
`*`(100) %>% # 百分比变化
sprintf("%.2f", .) %>%
as.data.frame() %>%
set_colnames(str_c("$\\theta=", sprintf("%.2f", theta), "$"))
}) %>%
mutate(
source_of_barrier = barrier_names,
parameters = barrier_parameters %>% sprintf("%.2f", .)
) %>%
select(source_of_barrier, parameters, everything())
table7 %>%
rename(
`Source of Barrier` = source_of_barrier,
`Estimated Geograpy Parameters` = parameters
) %>%
prettify(caption = "TABLE VII \\\n Geographic Barriers") %>%
footnote(
general = "The estimated parameters governing geographic barriers are the same as those shown in Table III. For an estimated parameter $\\hat{d}$, the implied effect on cost is $(e^{-\\hat{d}/\\theta}-1) \\times 100\\%$.",
general_title = "\n Note: ",
footnote_as_chunk = T
) %>%
add_header_above(c(
" " = 1, " " = 1,
"Implied Barrier's $\\%$ Effect on Cost" = 3
)) %>%
group_rows("Destination country:", 11, 29,
bold = FALSE,
label_row_css = "border-top: 1px solid;"
)| Source of Barrier | Estimated Geograpy Parameters | \(\theta=3.60\) | \(\theta=8.28\) | \(\theta=12.86\) |
|---|---|---|---|---|
| Distance [0, 375) | -3.10 | 136.51 | 45.39 | 27.25 |
| Distance [375, 750) | -3.66 | 176.74 | 55.67 | 32.97 |
| Distance [750, 1500) | -4.03 | 206.65 | 62.77 | 36.85 |
| Distance [1500, 3000) | -4.22 | 222.75 | 66.44 | 38.82 |
| Distance [3000, 6000) | -6.06 | 439.04 | 108.02 | 60.25 |
| Distance [6000, maximum) | -6.56 | 518.43 | 120.82 | 66.54 |
| Shared border | 0.30 | -7.89 | -3.51 | -2.27 |
| Shared language | 0.51 | -13.25 | -5.99 | -3.90 |
| Both in EEC/EU | 0.04 | -1.02 | -0.44 | -0.29 |
| Both in EFTA | 0.54 | -13.85 | -6.28 | -4.09 |
| Destination country: | ||||
| Australia | 0.24 | -6.35 | -2.81 | -1.82 |
| Austria | -1.68 | 59.37 | 22.46 | 13.94 |
| Belgium | 1.12 | -26.74 | -12.65 | -8.34 |
| Canada | 0.69 | -17.42 | -7.99 | -5.22 |
| Denmark | -0.51 | 15.15 | 6.33 | 4.03 |
| Finland | -1.33 | 44.88 | 17.49 | 10.94 |
| France | 0.22 | -5.90 | -2.61 | -1.69 |
| Germany | 1.00 | -24.27 | -11.39 | -7.49 |
| Greece | -2.36 | 92.45 | 32.93 | 20.11 |
| Italy | 0.07 | -1.97 | -0.86 | -0.56 |
| Japan | 1.59 | -35.62 | -17.43 | -11.60 |
| Netherlands | 1.00 | -24.33 | -11.42 | -7.51 |
| New Zealand | 0.07 | -1.83 | -0.80 | -0.52 |
| Norway | -1.00 | 32.06 | 12.85 | 8.10 |
| Portugal | -1.21 | 39.82 | 15.69 | 9.84 |
| Spain | -1.16 | 37.85 | 14.98 | 9.40 |
| Sweden | -0.02 | 0.69 | 0.30 | 0.19 |
| United Kingdom | 0.81 | -20.23 | -9.36 | -6.13 |
| United States | 2.46 | -49.49 | -25.70 | -17.40 |
|
Note: The estimated parameters governing geographic barriers are the same as those shown in Table III. For an estimated parameter \(\hat{d}\), the implied effect on cost is \((e^{-\hat{d}/\theta}-1) \times 100\%\). |
||||
estimation 部分，通过各种方式估计 \(\theta\)，符号一致，数值有大有小。
选择中间的、也是最简单的矩估计值 \(8.28\)，用于后面的 simulation