Using R for Trade Policy Analysis
R for Economic Research
Author: Nguyen Chi Dung
Motivation
Năm 2012, United Nations Conference On Trade And Development - UNCTAD và WTO phát hành Practical Guide to Trade Policy Analysis - một tài liệu kết hợp cả khía cạnh lí thuyết và phân tích thực chứng cho thương mại quốc tế. UNCTAD và WTO công khai đầy đủ Stata Codes (và cả dữ liệu) cho những kết quả được trình bày trong tập tài liệu này.
Sách này cùng với Stata Codes là một nguồn tài liệu tốt cho những ai theo đuổi bậc học sau đại học và tập trung vào việc ứng dụng những kiến thức và áp dụng của Kinh Tế Lượng (Econometrics) nhằm kiểm tra và đánh giá các lí thuyết kinh tế trong thương mại và mậu dịch quốc tế.
Sử dụng data được cung cấp và dựa trên Stata Codes được cung cấp, Massimiliano Porto viết cuốn Using R for Trade Policy Analysis tái hiện lại các phân tích và kết quả được trình bày trong tài liệu trên. Tuy nhiên textbook này có hai nhược điểm dưới đây:
Không nhất quán về coding và syntax (đọc thêm tại đây) gây khó khăn cho người học, đặc biệt là những người mới học R hoặc chưa có nhiều kinh nghiệm sử dụng R.
Một số R codes là sai. Chẳng hạn tính toán biến trễ của nhập khẩu (lag_totimppcode, page 21) là sai, dẫn đến cách kết quả khác cũng là sai.
Do vậy chuỗi bài giảng này được viết nhằm khắc phục hai nhược điểm trên. Cụ thể:
Phần xử lí dữ liệu sẽ được thực hiện nhất quán bằng gói dplyr còn phần hình ảnh hóa dữ liệu sẽ được thực hiện bằng gói ggplot2.
Sửa lại/hiệu chỉnh những đoạn R Codes mà tác giả Massimiliano Porto viết không đúng.
Chuỗi bài giảng này là phù hợp với những bạn muốn hoàn thiện kĩ năng xử lí số liệu, thực hiện các phân tích cho mục tiêu công bố nghiên cứu trên các tạp chí quốc tế bằng ngôn ngữ R. Cuốn sách này cũng phù hợp với những ai muốn chuyển dần các phân tích của mình từ Stata sang R.
Requirements
Người thực hành theo tài liệu này cần: (1) có kiến thức nền cơ bản về Kinh Tế Học cũng như Thống Kê/Kinh Tế Lượng, và (2) đã sử dụng R ở trình độ cơ bản.
Tự download tài liệu A Practical Guide to Trade Policy Analysis cùng dữ liệu (cũng như Stata Codes) tại đây.
Chapter 1: Analyzing Trade Flows
1. Overall openness
Trong mục này chúng ta sẽ tái lập lại Figure 1.1 nhằm đánh giá và ước lượng độ mở của nền kinh tế với bộ dữ liệu openness.dta. Mục tiêu thực hành là:
- Load dữ liệu từ file .dta (của phần mềm Stata).
- Tạo cột biến mới bằng hàm
dplyr::mutate(). - Chọn một số cột biến bằng hàm
dplyr::select().
Trước hết đọc bộ dữ liệu openness.dta bằng hàm read_dta() của gói haven như sau:
# Clear our workspace:
rm(list = ls())
# Load haven package:
library(haven)
# Load openness.dta:
openness <- read_dta("openness.dta")Xem qua dữ liệu:
# Some observations:
head(openness, 6)## # A tibble: 6 x 10
## ccode year openc openk pop gdp_current ldlock island remoteness
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AGO 1976 77.3 74.6 5943. NA 0 0 6524.
## 2 AGO 1977 78.6 74.7 6164. NA 0 0 4004.
## 3 AGO 1978 77.3 73.4 6287. NA 0 0 6895.
## 4 AGO 1979 78.3 72.4 6453. NA 0 0 4345.
## 5 AGO 1980 70.4 68.9 6743. NA 0 0 4385.
## 6 AGO 1981 84.9 80.7 6880. NA 0 0 5638.
## # ... with 1 more variable: remoteness_head <dbl>Xem một số thông tin về bộ dữ liệu này:
nrow(openness) # Number of rows or observation. ## [1] 3161ncol(openness) # Number of columns. ## [1] 10Dữ liệu ở dạng .dta thường có label mô tả sơ bộ ý nghĩa các cột biến. Chúng ta có thể xem qua mô tả này bằng hàm get_label() của thư viện sjlabelled như sau:
# Load sjlabelled package:
library(sjlabelled)
# Show labels:
get_label(openness)## ccode year
## "iso-code" "year"
## openc openk
## "trade openness in current terms" "trade openness in real terms"
## pop gdp_current
## "total population" "gdp in current terms"
## ldlock island
## "1 if landlock" "1 if island"
## remoteness remoteness_head
## "remoteness" "remoteness according to Head"Như vậy cột ldlock có mô tả là “1 if landlock”. Chúng ta có thể suy luận rằng đây là biến dummy ở đó 1 có nghĩa là quốc gia không có biển, và 0 tức là ngược lại.
Tạo một số biến số mới có kí hiệu và mô tả tại page 16 của tài liệu. Cụ thể:
- pop1000 (do đơn vị của pop hiện đang là 1000) bằng cách nhân 1000 với pop.
- gdppc (GDP per capita) bằng thương số của gdp_current chia cho pop1000.
- Logarit cơ số tự nhiên của gdppc, kí hiệu là ln_gdppc.
Sử dụng hàm mutate() của gói dplyr để tạo ra những biến này như sau:
# Load dplyr package:
library(dplyr)
# Create some new variables:
openness %>%
mutate(pop1000 = 1000*pop,
gdppc = gdp_current / pop1000,
ln_gdppc = log(gdppc)) -> opennessSử dụng hàm filter() để lấy ra chỉ các quan sát mà: (1) thuộc về năm 2000, và (2) có openc >= 200 như mô tả tại page 16:
openness %>% filter(year == 2000 & openc <= 200) -> openness_2000Bộ dữ liệu openness_2000 ở trên có thể được tạo ra bằng R codes sau (nếu bạn đã quen thuộc với ý tưởng của toán tử Pipe):
openness %>%
mutate(pop1000 = 1000*pop,
gdppc = gdp_current / pop1000,
ln_gdppc = log(gdppc)) %>%
filter(year == 2000 & openc <= 200) -> openness_2000Figure 1.1 tại page 16 (tương ứng với Stata Codes ở file openness.do) có thể được tạo ra bằng R codes như sau:
# Replicate Figure 1.1:
point_color <- "#149ccf" # Color code for point.
curve_color <- "#c099b8" # Color code for curve.
# Use Ubuntu Condensed font:
library(showtext)
showtext_auto()
# Select Ubuntu Condensed font for ploting:
font_add_google(name = "Ubuntu Condensed", family = "ubu")
my_font <- "ubu"
# Load ggplot2 package:
library(ggplot2)
# Left plot:
openness_2000 %>%
ggplot(aes(x = gdppc, y = openc, color = "Openness")) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, formula = y ~ x + I(x^2), aes(color = "Fitted values")) +
scale_color_manual(values = c(curve_color, point_color), name = "") +
theme_classic() +
theme(text = element_text(family = my_font)) +
theme(legend.position = "bottom") +
labs(x = "GDP per capita", y = "",
title = "Figure 1.1: Trade openness and GDP per capita, 2000",
subtitle = "(a) Quadratic fit") +
theme(plot.title = element_text(color = point_color, size = 16)) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
theme(axis.text = element_text(size = 12)) +
theme(legend.text = element_text(size = 10.5)) +
theme(legend.box.background = element_rect(colour = "black"),
legend.background = element_blank()) -> left_plot
# Right plot:
openness_2000 %>%
ggplot(aes(x = ln_gdppc, y = openc, color = "Openness")) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, formula = y ~ x + I(x^2), aes(color = "Fitted values")) +
scale_color_manual(values = c(curve_color, point_color), name = "") +
theme_classic() +
theme(text = element_text(family = my_font)) +
theme(legend.position = "bottom") +
labs(x = "Log GDP per capita", y = "",
title = "",
subtitle = "(b) Quadratic fit after log transformation") +
theme(plot.subtitle = element_text(hjust = 0.5)) +
theme(axis.text = element_text(size = 12)) +
theme(legend.text = element_text(size = 10.5)) +
theme(legend.box.background = element_rect(colour = "black"),
legend.background = element_blank()) -> right_plot
# Load gridExtra package:
library(gridExtra)
# Combine the two plots:
grid.arrange(left_plot, right_plot, ncol = 2)
2. Trade composition
Để tái tạo Figure 1.2 (page 19) trước hết chúng ta cần tính similarity index (page 20) như mô tả của Helpman (1987) từ bộ dữ liệu GravityData.dta. Dưới đây là R tính toán chỉ số SI này cho Đức ở năm 2004 (tương đương với Stata codes ở file overlap_trade.do):
# Load data (may be time-consuming):
haven::read_dta("GravityData.dta") -> gravity_df
# Remove duplications and calculate similarity index for German in 2004:
gravity_df %>%
filter(year == 2004 & ccode == "DEU") %>%
select(ccode, pcode, year, cgdp_c2000, pgdp_c2000) %>%
distinct() %>%
mutate(temp1 = cgdp_c2000 / (cgdp_c2000 + pgdp_c2000),
temp2 = pgdp_c2000 / (cgdp_c2000 + pgdp_c2000),
simil_index = 1 - temp1^2 - temp2^2)-> simil_for_DEU
# Remove missing data points:
simil_for_DEU %>% filter(!is.na(simil_index)) -> simil_for_DEUKế tiếp chúng ta tính toán Grubel-Lloyd (GL) index đo lường thương mại nội ngành (intra-industrytrade, page 19) của Đức từ bộ số liệu germany_trade_2004_hs6.dta:
#====================================================
# Open 2004 (hs6) trade data of Germany,
# reshape and modify it for future use
#====================================================
# Load data:
haven::read_dta("germany_trade_2004_hs6.dta") -> german_trade_hs6
# Load stringr package for text processing:
library(stringr)
german_trade_hs6 %>%
rename(ccode = reporter, pcode = partner) %>% # Rename for some columns.
mutate(flow_name = case_when(flow_name == "Gross Exp." ~ "Exports", TRUE ~ "Imports")) %>%
mutate(product = str_squish(product)) %>%
mutate(n_text = str_count(product)) %>%
filter(n_text >= 6) %>%
mutate(n_text = NULL, nomen = NULL, rownum = NULL) -> germman_trade_long
# Convert to wide form:
library(tidyr)
germman_trade_long %>% tidyr::spread(key = "flow_name", value = "trade_value") -> overlap_temp
# Replace NA by zero and calculate Grubel-Lloyd (GL) index by product:
overlap_temp %>%
mutate(Exports = replace_na(Exports, 0),
Imports = replace_na(Imports, 0)) -> overlap_temp
overlap_temp %>%
group_by(pcode, product) %>%
summarise(gl_i_j_k = 1 - abs((Exports - Imports)) / (Exports + Imports)) %>%
mutate(gl_i_j_k = abs(gl_i_j_k)) %>%
ungroup() -> df_gl
# Join the two data sets:
full_join(overlap_temp, df_gl, by = c("pcode", "product")) -> overlap_tempTính trade overlap index (page 20):
overlap_temp %>%
group_by(pcode) %>%
summarise(x1 = sum(Exports), x2 = sum(Imports)) %>%
mutate(denom = x1 + x2) -> d1
overlap_temp %>%
filter(gl_i_j_k > 0) %>%
group_by(pcode) %>%
summarise(x11 = sum(Exports), x22 = sum(Imports)) %>%
mutate(numer = x11 + x22) -> d2
full_join(d1, d2, by = c("pcode")) %>%
mutate(overlap = numer / denom) -> d3
full_join(d3 %>% select(pcode, overlap), simil_for_DEU, by = c("pcode")) -> overlap_data_DEUR codes cho Fugure 1.2:
# Figure 1.2:
library(ggrepel)
overlap_data_DEU %>%
ggplot(aes(simil_index, overlap)) +
geom_point(aes(color = "Openness")) +
geom_smooth(method = "lm", se = FALSE, formula = y ~ x + I(x^2), aes(color = "Fitted values"), size = 1.2) +
geom_text_repel(aes(label = pcode), size = 3.5, family = my_font) +
scale_color_manual(values = c(curve_color, point_color), name = "") +
theme_classic() +
theme(text = element_text(family = my_font)) +
theme(legend.position = "bottom") +
labs(x = "Similarity index", y = "",
title = "Figure 1.2: Overlap trade and country-similarity index vis-à-vis germany, 2004",
caption = "Author calculations from World Bank WDI and UN Comtrade") +
theme(plot.title = element_text(color = point_color, size = 16)) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
theme(axis.text = element_text(size = 12)) +
theme(legend.text = element_text(size = 10.5)) +
theme(legend.box.background = element_rect(colour = "black"),
legend.background = element_blank()) 
(Đang được viết tiếp…)