-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path4.0 analysis_medium_functions.R
202 lines (175 loc) · 9.03 KB
/
4.0 analysis_medium_functions.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
#' Analyse medium-term forecasts bias
#'
#' @param data df with forecasts at different horizons and actual values (named respectively variable7/8/9/10/11/12
#' and targety_first)
#' @param regressions character vector. Different regressions to perform.
#' @param output_type character string. Can be "appendix_table", "share_plot","magnitude_table",
#' "share_plot_geo" or "magnitude_table_geo".
#' @param export_path character string. Path to export table or graph.
#' @return Output depends on parameter output type.
#' The first produces a tex table with value of the intercept for each country and forecast horizon.
#' The second produces two plots with the share of countries with significant intercept
#' at each horizon (dividing between Spring and Fall.)
#' The third produces a table with the summary statistics of biases for each horizon and type
#' of bias i.e optimistic or pessimistic. For each horizon, Spring and Fall issues are pooled together.
#' The fourth produces the same plots as the second, but dividing by WEO geographical group.
#' The fifth produces the same table as the third, but dividing by WEO geographical group.
#'
analyse_medium_bias <- function(data, regressions, output_type, export_path){
# Run regressions for each country and formula, produces nested list:
regressions <- data %>%
mutate_at(vars(contains("variable")),funs(targety_first - .)) %>%
split(.$country_code) %>%
map( ~ map(regressions, function(x){
tryCatch(lm(x, .x), error = function(e){
cat(crayon::red("Could not run the regression. Check data\n"))
})
}))
# Wrangle in dataframe format with estimate, country and horizon:
df_bias <- regressions %>%
modify_depth(2, ~ tryCatch(summary(.x), error = function(e){
cat(crayon::red("Could not run regression. Check dataframe \n"))
})) %>%
modify_depth(2, ~ tryCatch(.x[["coefficients"]], error = function(e){
cat(crayon::red("Could not run regression. Check dataframe \n"))
})) %>%
map(~ discard(.x, ~ length(.x) != 4)) %>%
modify_depth(2, ~ as.data.frame(.x)) %>%
map(~ bind_rows(.x, .id = "horizon")) %>%
bind_rows(.id = "country_code") %>%
mutate(issue = case_when(as.numeric(horizon) %% 2 == 1 ~ "Fall",
T ~ "Spring"),
horizon = case_when(horizon == 1 | horizon == 2 ~ "H=3",
horizon == 3 | horizon == 4 ~ "H=4",
horizon == 5 | horizon == 6 ~ "H=5")) %>%
mutate(Estimate = round(Estimate,2)) %>%
mutate(Estimate = case_when(`t value` > 1.96 | `t value` < -1.96 ~ str_replace(as.character(Estimate), "$","**"),
(`t value` > 1.68 & `t value` < 1.96) | (`t value` < -1.68 & `t value` > -1.96) ~ str_replace(as.character(Estimate), "$", "*"),
TRUE ~ as.character(Estimate))) %>%
mutate(country = countrycode(country_code,"imf","country.name")) %>%
select(country_code, country, horizon, issue, Estimate)
# Produces output:
if(output_type == "appendix_table"){
df_bias %>%
unite("horizon",horizon, issue, sep = ",") %>%
select(-country_code) %>%
spread(horizon, Estimate) %>%
rename(Country = country) %>%
stargazer(summary = F,
rownames = F,
out = export_path)
}
else if(output_type == "share_plot"){
share_aggregate <- df_bias %>%
split(.$horizon) %>%
map(~ .x %>% mutate(negative_significant = case_when(str_detect(Estimate, "\\*") & str_detect(Estimate, "-") ~ 1,
T ~ 0))) %>%
map(~ .x %>% mutate(positive_significant = case_when(str_detect(Estimate, "\\*") & str_detect(Estimate, "-", negate = T) ~ 1,
T ~ 0))) %>%
map(~ .x %>% split(.$issue)) %>%
modify_depth(2, ~ .x %>% summarise_at(vars(contains("significant")), mean, na.rm = T)) %>%
map(~ .x %>% bind_rows(.id = "issue")) %>%
bind_rows(.id = "horizon") %>%
gather("sign","share",negative_significant:positive_significant) %>%
mutate(sign = case_when(sign == "negative_significant" ~ "Optimistic",
T ~ "Pessimistic")) %>%
split(.$issue) %>%
map(~ .x %>%
ggplot(aes(sign, share, fill = sign)) +
geom_col(width = 0.4) +
geom_text(aes(label = round(share,2)), size = 5, vjust = -0.5) +
facet_wrap(~ horizon) +
theme_minimal() +
ylim(0,1) +
xlab("") +
ylab("Share of countries (%)") +
scale_fill_grey() +
theme(legend.position = "bottom") +
theme(strip.text.x = element_text(size = 16),
axis.text.x = element_blank(),
axis.text.y = element_text(size = 18),
axis.title = element_text(size = 21),
legend.text = element_text(size = 16)) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()) +
labs(fill = ""))
# Export:
share_aggregate %>%
iwalk(~ ggsave(filename = paste0(export_path,.y,".pdf"),.x))
}
else if(output_type == "magnitude_table"){
table_magnitude <- df_bias %>%
split(.$horizon) %>%
map(~ .x %>% filter(str_detect(Estimate,"\\*"))) %>%
map(~ .x %>% mutate(negative = case_when(str_detect(Estimate,"-") ~ 1,
T ~ 0))) %>%
map(~ .x %>% mutate(Estimate = as.numeric(str_remove(Estimate, "\\*+")))) %>%
map(~ .x %>% group_by(negative) %>% summarise(mean_bias = round(mean(Estimate, na.rm = T),2),
median_bias = round(median(Estimate, na.rm = T),2),
max_bias = round(max(Estimate),2),
min_bias = round(min(Estimate),2))) %>%
bind_rows(.id = "horizon") %>%
mutate(negative = case_when(negative == 1 ~ "Optimistic",
T ~ "Pessimistic")) %>%
arrange(negative) %>%
setNames(c("Horizon","Type of bias","Mean","Median", "Min.", "Max."))
table_magnitude %>%
stargazer(summary = F,
rownames = F,
out = export_path)
}
else if(output_type == "share_plot_geo"){
share_aggregate_group <- df_bias %>%
merge(geo_group,by=c("country_code")) %>%
split(.$horizon) %>%
map(~ .x %>% mutate(negative_significant = case_when(str_detect(Estimate, "\\*") & str_detect(Estimate, "-") ~ 1,
T ~ 0))) %>%
map(~ .x %>% mutate(positive_significant = case_when(str_detect(Estimate, "\\*") & str_detect(Estimate, "-", negate = T) ~ 1,
T ~ 0))) %>%
map(~ .x %>% split(.$issue)) %>%
modify_depth(2, ~ .x %>% group_by(group)) %>%
modify_depth(2, ~ .x %>% summarise_at(vars(contains("significant")), mean, na.rm = T)) %>%
map(~ .x %>% bind_rows(.id = "issue")) %>%
bind_rows(.id = "horizon") %>%
gather("sign","share",negative_significant:positive_significant) %>%
mutate(sign = case_when(sign == "negative_significant" ~ "Optimistic",
T ~ "Pessimistic")) %>%
split(.$issue) %>%
map(~ .x %>%
ggplot(aes(group,share, fill = sign)) +
geom_col(position = "dodge",width = 0.5) +
coord_flip() +
facet_wrap(~horizon) +
theme_minimal() +
xlab("") +
ylab("Share of countries (%)") +
scale_fill_grey() +
theme(legend.position = "bottom") +
labs(fill="") +
theme(strip.text.x = element_text(size = 16),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 18),
axis.title = element_text(size = 21),
legend.text = element_text(size = 16))
)
share_aggregate_group %>%
iwalk(~ ggsave(filename = paste0(export_path,.y,"_group.pdf"),.x))
}
else if(output_type == "magnitude_table_geo"){
df_bias %>%
merge(geo_group,by=c("country_code")) %>%
split(.$horizon) %>%
map(~ .x %>% filter(str_detect(Estimate,"\\*"))) %>%
map(~ .x %>% filter(str_detect(Estimate,"-"))) %>%
map(~ .x %>% group_by(group)) %>%
map(~ .x %>% mutate(Estimate = as.numeric(str_remove_all(Estimate,"\\*")))) %>%
map(~ .x %>% summarise(mean = round(mean(Estimate, na.rm = T),2),
median = round(median(Estimate, na.rm = T),2),
max = round(min(Estimate, na.rm = T),2))) %>%
bind_rows(.id = "horizon") %>%
setNames(c("Horizon","Geo. group","Mean","Median","Max.")) %>%
stargazer(summary = F,
rownames = F,
out = export_path)
}
}