Appendix A: Further Analysis on the Trade Networks
First we will be looking at the results of the hierarchical cluster, through the use of dendrograms.
Figure 27: Dendrogram of the Import Network
From figure 27, we can see the various groupings of the import network through the use of hierarchical cluster. If we cut the graph at the height of 1.1, we have four groups of countries.
The first group of countries goes from Austria to Poland. The second group of countries goes from Belarus to Bulgaria. The third group of countries goes from Canada to Slovenia. The fourth group of countries goes from Turkey to Malta. The geographical effect is still somehow observable in group one, e.g. Norway and Sweden are grouped close together, but the geographical effect is not as pronounced in the successive groups. It is interesting to note that Germany and the Netherlands are grouped together.
It is interesting to note that group three and four mainly consist of import partners for the EU, meaning that these countries supply the paper and cardboard waste for the EU.
Austria Slovakia Hungary Italy Netherlands Germany Estonia Russia Lithuania Latvia Sweden Norway Finland Poland Belarus Moldova Montenegro Bosnia Serbia Romania Macedonia Croatia Ukraine Czechia Albania Bulgaria Canada Dominican Rep Costa Rica France Spain Belgium Denmark Slovenia Turkey Andorra Gibraltar Morocco Melilla Algeria United States Portugal Switzerland Cyprus Israel Trinidad Thailand Ireland Iceland Greece Luxembourg Malta
0.00.20.40.60.81.01.2
Dendrogram of the Import Network
hclust (*, "complete") Source: Eurostat
Height
Figure 28: Dendrogram of the Export Network
If we cut the graph at the height of 1.2 in figure 28, two distinctive groups are observed. The first group goes from Hungary until Panama, and the second group goes from Poland to South Korea. The geographical effect is still there, but only at the lower levels of clustering. It is interesting to note that Germany and the Netherlands are still grouped together.
Hungary Slovakia Czechia Croatia Slovenia Austria Romania Bulgaria Serbia Chile Belarus Ukraine Estonia Lithuania Russia Israel Lebanon Saudi Arabia Bangladesh Moldova Macedonia Ecuador Faroe Islands Norway Colombia Mexico Cuba Morocco Tunisia Djibouti Hong Kong Egypt Ghana Panama Poland Latvia Finland Denmark Sweden Syria Netherlands Germany Italy France Portugal Belgium Cyprus Greece Ireland Malta Australia China Switzerland Spain Luxembourg Philippines Japan Turkey India Pakistan Vietnam Taiwan Malaysia Indonesia Thailand United States Singapore South Korea
0.00.20.40.60.81.01.21.4
Dendrogram of the Export Network
hclust (*, "complete") Source: Eurostat
Height
Figure 29: Dendrogram of the Full Trade Network
Figure 29 shows the clustering result of the full trade network of paper and cardboard waste.
Compared to the clustering results of the export network, the results of the full network remains largely unchanged. This is unsurprising as the clustering algorithm takes into account of the quantity traded, and the quantity exported is much larger than the quantity imported thus resulting in two similar graphs.
Hungary Slovakia Czechia Croatia Slovenia Austria Romania Bulgaria Serbia Chile Belarus Ukraine Estonia Lithuania Russia Israel Lebanon Saudi Arabia Bangladesh Moldova Macedonia Ecuador Faroe Islands Norway Colombia Mexico Cuba Morocco Tunisia Djibouti Hong Kong Egypt Ghana Panama Poland Latvia Finland Denmark Sweden Syria Netherlands Germany Italy France Portugal Belgium Cyprus Greece Ireland Malta Australia China Switzerland Spain Luxembourg Philippines Japan Turkey India Pakistan Vietnam Taiwan Malaysia Indonesia Thailand United States Singapore South Korea
0.00.20.40.60.81.01.21.4
Dendrogram of the Full Trade Network
hclust (*, "complete") Source: Eurostat
Height
Appendix B: R Codes Used in this Thesis
# Thesis Chapter 1 1 # Municipal Waste 2
3
require(readxl) 4
require(tidyverse) 5
require(ggpubr) 6
require(viridis) 7 require(plm) 8
require(lmtest) 9
require(tseries) 10
11
#### International Comparison ####
12 13
# Total Municipal Waste 14
total_desc <- read_excel("Desktop/Chapter 1/descriptive/total desc.xlsx") 15
totall = total_desc %>% rename(Region = country) %>% pivot_longer(-c(Region), names_to 16 = "Year", values_to = "total")
17 18
ggplot(data = totall, aes(x=Year,y=total,group=Region)) + 19
geom_point(aes(color=Region,shape=Region),size=4) + 20
labs(x="Year",y="Total Municipal Waste in kg per capita", title = "Total Municipal Waste 21 Generated over Time", caption="Source: Eurostat") +
22
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 23
24
# Landfilled Municipal Waste 25
landfill_desc <- read_excel("Desktop/Chapter 1/descriptive/landfill desc.xlsx") 26
landfilledl = landfill_desc %>% rename(Region = country) %>% pivot_longer(-c(Region), 27 names_to = "Year", values_to = "landfill")
28 29
ggplot(data = landfilledl, aes(x=Year,y=landfill,group=Region)) + 30 geom_point(aes(color=Region,shape=Region),size=4) +
31
labs(x="Year",y="Landfilled Municipal Waste in kg per capita", title = "Landfilled 32
Municipal Waste over Time", caption="Source: Eurostat") + 33
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 34
35
# Incinerated Municipal Waste
36 incinerated_desc <- read_excel("Desktop/Chapter 1/descriptive/incinerated desc.xlsx") 37
incineratedl = incinerated_desc %>% rename(Region = country) %>% pivot_longer(- 38
c(Region), names_to = "Year", values_to = "incinerate") 39
40
ggplot(data = incineratedl, aes(x=Year,y=incinerate,group=Region)) + 41 geom_point(aes(color=Region,shape=Region),size=4) +
42
labs(x="Year",y="Incinerated Municipal Waste in kg per capita", title = "Incinerated 43
Municipal Waste for Disposal over Time", caption="Source: Eurostat") + 44
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 45
46 # Recovered Municipal Waste 47
energy_recovery_desc <- read_excel("Desktop/Chapter 1/descriptive/energy recovery 48
desc.xlsx") 49
recoveredl = energy_recovery_desc %>% rename(Region = country) %>% pivot_longer(- 50 c(Region), names_to = "Year", values_to = "recover")
51
52
ggplot(data = recoveredl, aes(x=Year,y=recover,group=Region)) + 53
geom_point(aes(color=Region,shape=Region),size=4) +
54 labs(x="Year",y="Recovered Municipal Waste in kg per capita", title = "Recovered 55
Municipal Waste for Energy over Time", caption="Source: Eurostat") + 56
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 57
58 # Recycled Municipal Waste 59
recyc_desc <- read_excel("Desktop/Chapter 1/descriptive/recyc desc.xlsx") 60
recyclel = recyc_desc %>% rename(Region = country) %>% pivot_longer(-c(Region), 61
names_to = "Year", values_to = "recycle") 62
63 ggplot(data = recyclel, aes(x=Year,y=recycle,group=Region)) + 64
geom_point(aes(color=Region,shape=Region),size=4) + 65
labs(x="Year",y="Recycled Municipal Waste in kg per capita", title = "Recycled Municipal 66
Waste of Materials over Time", caption="Source: Eurostat") + 67
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 68
69 rm(list = ls()) 70
#### Preparing Data for Panel Regression ####
71 72
income <- read_excel("Desktop/Chapter 1/FCE Households.xlsx")
73 incomelong = income %>% pivot_longer(-c(country), names_to = "year", values_to = "fceo") 74
75
total <- read_excel("Desktop/Chapter 1/total.xlsx") 76
totallong = total %>% pivot_longer(-c(country), names_to = "year", values_to = "totalo") 77
totalw = left_join(totallong, incomelong, by = c("country","year")) 78 saveRDS(totalw, file = "totalw.rds")
79 80
landfilled <- read_excel("Desktop/Chapter 1/landfilled.xlsx") 81
landfilledlong = landfilled %>% pivot_longer(-c(country), names_to = "year", values_to = 82
"lando")
83 landw = left_join(landfilledlong, incomelong, by = c("country","year")) 84
saveRDS(landw, file = "landw.rds") 85
86
recycling_of_materials <- read_excel("Desktop/Chapter 1/recycling of materials.xlsx") 87
recyclong = recycling_of_materials %>% pivot_longer(-c(country), names_to = "year", 88 values_to = "recyco")
89
recycw = left_join(recyclong, incomelong, by = c("country","year")) 90
saveRDS(recycw, file = "recycw.rds") 91
92 #### Panel Descriptives ####
93 94
# TOTAL 95
totalw <- readRDS("~/Desktop/Chapter 1/totalw.rds") 96
summary(totalw) 97 98
panel1 = totalw %>% mutate(fce = log(fceo),fce2 = 99
log(fceo)*log(fceo),fce3=log(fceo)*log(fceo)*log(fceo),total=log(totalo),.keep="unused") 100
101
ggplot(data = panel1, aes(x=fce,y=total)) + 102
geom_point(aes(color=country)) + 103
labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 104
PPS per capita",y="ln(Total Municipal Waste) in kg per capita", title = "Total Municipal 105
Waste vs Income", caption="Source: Eurostat") + theme(legend.position = "bottom") + 106 scale_color_viridis(discrete=TRUE)
107 108
ggplot(data = panel1, aes(x=fce,y=total)) + 109
geom_point(aes(color=country)) +
110 labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 111
PPS per capita",y="ln(Total Municipal Waste) in kg per capita", title = "Total Municipal 112
Waste vs Income by Country", caption="Source: Eurostat") + theme(legend.position = 113
"bottom") + scale_color_viridis(discrete=TRUE) + 114
facet_wrap(~country,scales = "free") + 115 theme(legend.position = "none")
116 117
# LANDFILL 118
119
landw <- readRDS("~/Desktop/Chapter 1/landw.rds") 120
summary(landw) 121 122
panel2 = landw %>% mutate(fce = log(fceo),fce2 = 123
log(fceo)*log(fceo),fce3=log(fceo)*log(fceo)*log(fceo),land=log(lando),.keep="unused") 124
%>% mutate_if(is.numeric, ~ifelse(abs(.) == Inf,0,.)) 125 126
ggplot(data = panel2, aes(x=fce,y=land)) + 127
geom_point(aes(color=country)) + 128
labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 129
PPS per capita",y="ln(Landfilled Municipal Waste) in kg per capita", title = "Landfilled 130 Municipal Waste vs Income", caption="Source: Eurostat") + theme(legend.position = "bottom") 131
+ scale_color_viridis(discrete=TRUE) 132
133
ggplot(data = panel2, aes(x=fce,y=land)) + 134
geom_point(aes(color=country)) +
135 labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 136
PPS per capita",y="ln(Landfilled Municipal Waste) in kg per capita", title = "Landfilled 137
Municipal Waste vs Income by Country", caption="Source: Eurostat") + theme(legend.position 138
= "bottom") + scale_color_viridis(discrete=TRUE) + 139
facet_wrap(~country,scales = "free") + 140 theme(legend.position = "none")
141 142
# RECYCLING 143
144 recycw <- readRDS("~/Desktop/Chapter 1/recycw.rds") 145
summary(recycw) 146
147
panel3 = recycw %>% mutate(fce = log(fceo),fce2 = 148
log(fceo)*log(fceo),fce3=log(fceo)*log(fceo)*log(fceo),recyc=log(recyco),.keep="unused") 149 %>% mutate_if(is.numeric, ~ifelse(abs(.) == Inf,0,.))
150 151
ggplot(data = panel3, aes(x=fce,y=recyc)) + 152
geom_point(aes(color=country)) + 153
labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 154
PPS per capita",y="ln(Recycled Municipal Waste) in kg per capita", title = "Recycled 155
Materials of Municipal Waste vs Income", caption="Source: Eurostat") + theme(legend.position 156
= "bottom") + scale_color_viridis(discrete=TRUE) 157
158 ggplot(data = panel3, aes(x=fce,y=recyc)) + 159
geom_point(aes(color=country)) + 160
labs(color="Country", x="ln(Final Consumption Expenditure of Households) in current prices 161
PPS per capita",y="ln(Recycled Municipal Waste) in kg per capita", title = "Recycled 162 Materials of Municipal Waste vs Income by Country", caption="Source: Eurostat") + 163
theme(legend.position = "bottom") + scale_color_viridis(discrete=TRUE) + 164
facet_wrap(~country,scales = "free") + 165
theme(legend.position = "none") 166
167 #### TOTAL Panel Regression PANEL 1 ####
168 169
plm::is.pbalanced(panel1$country,panel1$year) # check if it is a balanced panel 170
171
# fixed effects 172
fe1 = plm(total ~ fce, data = panel1,index=c("country","year"),model="within") # total 173 summary(fe1, vcovHC)
174 175
# fixed effects with squared term 176
fe2 = plm(total ~ fce + fce2, data = panel1,index=c("country","year"),model="within") 177 summary(fe2, vcovHC)
178 179
# fixed effects with cubic term 180
fe3 = plm(total ~ fce + fce2 + fce3, data = 181
panel1,index=c("country","year"),model="within") 182 summary(fe3, vcovHC)
183 184
# random effects 185
re1 = plm(total ~ fce, data = panel1,index=c("country","year"),model="random") 186
summary(re1, vcovHC) 187 188
#random effects with squared term 189
re2 = plm(total ~ fce + fce2, data = panel1,index=c("country","year"),model="random") 190
summary(re2, vcovHC) 191
192 #random effects with cubic term 193
re3 = plm(total ~ fce + fce2 + fce3, data = 194
panel1,index=c("country","year"),model="random") 195
summary(re3, vcovHC) 196 197
# fixed effects with time 198
fet1 = plm(total ~ fce + factor(year), data=panel1, model="within") 199
summary(fet1,vcovHC) 200
201 fet2 = plm(total ~ fce + fce2 + factor(year), data=panel1, model="within") 202
summary(fet2,vcovHC) # with squared term.
203 204
fet3 = plm(total ~ fce + fce2 + fce3 + factor(year), data=panel1, model="within") 205
summary(fet3,vcovHC) # with cubic term 206
207 # random effects with time 208
ret1 = plm(total ~ fce + factor(year), data=panel1, 209
model="random",index=c("country","year")) 210
summary(ret1,vcovHC) 211 212
ret2 = plm(total ~ fce + fce2 + factor(year), data=panel1, 213
model="random",index=c("country","year")) 214
summary(ret2,vcovHC) # with squared term 215 216
ret3 = plm(total ~ fce + fce2 + fce3 + factor(year), data=panel1, 217
model="random",index=c("country","year")) 218
summary(ret3,vcovHC) # with cubic term 219
220 # panel diagnostics 221
222
cor(panel1$fce,panel1$total, method="pearson") # The correlation between the two 223
variables is 0.5564985.
224 225
pt1 = plm(total ~ fce + fce2 + fce3 + factor(year) + factor(country), data=panel1, 226 model="pooling")
227
plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 228
pooled OLS is not appropriate 229
pt2 = plm(total ~ fce + fce2 + factor(year) + factor(country), data=panel1, 230 model="pooling")
231
plmtest(pt2, type=c("bp")) # there is a panel effect 232
pt3 = plm(total ~ fce + factor(year) + factor(country), data=panel1, model="pooling") 233
plmtest(pt3, type=c("bp")) # there is a panel effect 234
235 pFtest(fet2,fe2) # reject HO, and state there is a time effect 236
pFtest(ret2,re2) # there is a time effect 237
plmtest(fet2, c("time"), type=("bp")) # there is a time effect 238
plmtest(ret2, c("time"), type=("bp")) # there is a time effect 239
240 adf.test(panel1$total, k=2) 241
# Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 242
series is stationary 243
adf.test(panel1$fce, k=2) # stationary 244
adf.test(panel1$fce2, k=2) # stationary 245 adf.test(panel1$fce3, k=2) # stationary 246
247
# Model diagnostics 248
# All of these tests need to be performed individually for the models we would like to 249 compare
250 251
phtest(fet2,ret2) # fail to reject null (random) and state that the random effects is more 252
appropriate 253
254 pcdtest(fet2) # not cross-sectional dependent 255
pcdtest(ret2) # not cross-sectional dependent 256
257
pbgtest(fet2) # there is serial correlation 258
pbgtest(ret2) # there is serial correlation 259
260
bptest(total ~ fce + factor(year), data = panel1, studentize=F) # fail to reject HO, thus 261
there's no heteroskedaticity 262
bptest(total ~ fce + fce2 + fce3 + factor(year), data = panel1, studentize=F) # reject HO, 263 thus there is heteroskedaticity
264 265
#### LANDFILLED Panel Regression PANEL 2 ####
266
267 plm::is.pbalanced(panel2$country,panel2$year) # check if it is a balanced panel 268
269
# fixed effects 270
fe1 = plm(land ~ fce, data = panel2,index=c("country","year"),model="within") # land 271
summary(fe1, vcovHC) 272 273
# fixed effects with squared term 274
fe2 = plm(land ~ fce + fce2, data = panel2,index=c("country","year"),model="within") 275
summary(fe2, vcovHC) 276
277
# fixed effects with cubic term
278 fe3 = plm(land ~ fce + fce2 + fce3, data = 279
panel2,index=c("country","year"),model="within") 280
summary(fe3, vcovHC) 281
282 # random effects 283
re1 = plm(land ~ fce, data = panel2,index=c("country","year"),model="random") 284
summary(re1, vcovHC) 285
286
#random effects with squared term
287 re2 = plm(land ~ fce + fce2, data = panel2,index=c("country","year"),model="random") 288
summary(re2, vcovHC) 289
290
#random effects with cubic term 291
re3 = plm(land ~ fce + fce2 + fce3, data = 292 panel2,index=c("country","year"),model="random") 293
summary(re3, vcovHC) 294
295
# fixed effects with time 296
fet1 = plm(land ~ fce + factor(year), data=panel2, model="within") 297 summary(fet1,vcovHC)
298 299
fet2 = plm(land ~ fce + fce2 + factor(year), data=panel2, model="within") 300
summary(fet2,vcovHC) # with squared term.
301 302
fet3 = plm(land ~ fce + fce2 + fce3 + factor(year), data=panel2, model="within") 303
summary(fet3,vcovHC) # with cubic term 304
305
# random effects with time
306 ret1 = plm(land ~ fce + factor(year), data=panel2, 307
model="random",index=c("country","year")) 308
summary(ret1,vcovHC) 309
310
ret2 = plm(land ~ fce + fce2 + factor(year), data=panel2, 311
model="random",index=c("country","year")) 312 summary(ret2,vcovHC) # with squared term 313
314
ret3 = plm(land ~ fce + fce2 + fce3 + factor(year), data=panel2, 315
model="random",index=c("country","year")) 316 summary(ret3,vcovHC) # with cubic term 317
318
# panel diagnostics 319
320 cor(panel2$fce,panel2$land, method="pearson") # The correlation between the two 321
variables is 0.5564985.
322 323
pt1 = plm(land ~ fce + fce2 + fce3 + factor(year) + factor(country), data=panel2, 324
model="pooling")
325 plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 326
pooled OLS is not appropriate 327
pt2 = plm(land ~ fce + fce2 + factor(year) + factor(country), data=panel2, 328
model="pooling") 329
plmtest(pt2, type=c("bp")) # there is a panel effect 330
pt3 = plm(land ~ fce + factor(year) + factor(country), data=panel2, model="pooling") 331 plmtest(pt3, type=c("bp")) # there is a panel effect
332 333
pFtest(fet2,fe2) # reject HO, and state there is a time effect 334
pFtest(ret2,re2) # there is a time effect 335 pFtest(fet3,fe3) # there is a time effect 336
pFtest(fet1,fe1) # there is a time effect 337
pFtest(ret3,re3) # there is a time effect 338
plmtest(fet2, c("time"), type=("bp")) # there is a time effect 339
plmtest(ret2, c("time"), type=("bp")) # there is a time effect 340 plmtest(fet3, c("time"), type=("bp")) # there is a time effect 341
342
adf.test(panel2$land, k=2) 343
# Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 344
series is stationary
345 adf.test(panel2$fce, k=2) # stationary 346
adf.test(panel2$fce2, k=2) # stationary 347
adf.test(panel2$fce3, k=2) # stationary 348
349
# Model diagnostics
350 # All of these tests need to be performed individually for the models we would like to 351
compare 352
353
phtest(fet2,ret2) # fail to reject null (random) and state that the random effects are more 354 appropriate
355
phtest(fet3,ret3) # random effects are more appropriate 356
357
pcdtest(fet2) # not cross-sectional dependent 358
pcdtest(ret2) # not cross-sectional dependent 359 pcdtest(fet3) # not cross-sectional dependent 360
pcdtest(ret3) # not cross-sectional dependent 361
362
pbgtest(fet2) # there is serial correlation 363
pbgtest(fet3) # there is serial correlation 364
pbgtest(ret2) # there is serial correlation 365 pbgtest(ret3) # there is serial correlation 366
367
bptest(land ~ fce + factor(year), data = panel2, studentize=F) # there is heteroskedaticity 368
bptest(land ~ fce + fce2 + factor(year), data = panel2, studentize=F) # reject HO, thus 369 there is heteroskedaticity
370
bptest(land ~ fce + fce2 + fce3 + factor(year), data = panel2, studentize=F) # reject HO, 371
thus there is heteroskedaticity 372
373 bptest(ret3, studentize=F) 374
375
#### Recycled Panel Regression PANEL 3 ####
376 377
plm::is.pbalanced(panel3$country,panel3$year) # check if it is a balanced panel 378 379
# fixed effects 380
fe1 = plm(recyc ~ fce, data = panel3,index=c("country","year"),model="within") # recyc 381
summary(fe1, vcovHC) 382
383
# fixed effects with squared term
384 fe2 = plm(recyc ~ fce + fce2, data = panel3,index=c("country","year"),model="within") 385
summary(fe2, vcovHC) 386
387
# fixed effects with cubic term
388 fe3 = plm(recyc ~ fce + fce2 + fce3, data = 389
panel3,index=c("country","year"),model="within") 390
summary(fe3, vcovHC) 391
392
# random effects
393 re1 = plm(recyc ~ fce, data = panel3,index=c("country","year"),model="random") 394
summary(re1, vcovHC) 395
396
#random effects with squared term 397
re2 = plm(recyc ~ fce + fce2, data = panel3,index=c("country","year"),model="random") 398 summary(re2, vcovHC)
399 400
#random effects with cubic term 401
re3 = plm(recyc ~ fce + fce2 + fce3, data = 402
panel3,index=c("country","year"),model="random") 403 summary(re3, vcovHC)
404 405
# fixed effects with time 406
fet1 = plm(recyc ~ fce + factor(year), data=panel3, model="within") 407 summary(fet1,vcovHC)
408 409
fet2 = plm(recyc ~ fce + fce2 + factor(year), data=panel3, model="within") 410
summary(fet2,vcovHC) # with squared term.
411
412 fet3 = plm(recyc ~ fce + fce2 + fce3 + factor(year), data=panel3, model="within") 413
summary(fet3,vcovHC) # with cubic term 414
415
# random effects with time 416
ret1 = plm(recyc ~ fce + factor(year), data=panel3, 417
model="random",index=c("country","year")) 418 summary(ret1,vcovHC)
419
420
ret2 = plm(recyc ~ fce + fce2 + factor(year), data=panel3, 421
model="random",index=c("country","year")) 422 summary(ret2,vcovHC) # with squared term 423
424
ret3 = plm(recyc ~ fce + fce2 + fce3 + factor(year), data=panel3, 425
model="random",index=c("country","year")) 426 summary(ret3,vcovHC) # with cubic term 427
428
# panel diagnostics 429
430
cor(panel3$fce,panel3$recyc, method="pearson") # The correlation between the two 431 variables is 0.5564985.
432 433
pt1 = plm(recyc ~ fce + fce2 + fce3 + factor(year) + factor(country), data=panel3, 434
model="pooling") 435
plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 436
pooled OLS is not appropriate
437 pt2 = plm(recyc ~ fce + fce2 + factor(year) + factor(country), data=panel3, 438
model="pooling") 439
plmtest(pt2, type=c("bp")) # there is a panel effect 440
pt3 = plm(recyc ~ fce + factor(year) + factor(country), data=panel3, model="pooling") 441 plmtest(pt3, type=c("bp")) # there is a panel effect
442 443
pFtest(ret1,re1) # there is a time effect 444
pFtest(fet3,fe3) # there is a time effect 445
pFtest(fet1,fe1) # there is a time effect 446 pFtest(ret3,re3) # there is a time effect 447
plmtest(fet1, c("time"), type=("bp")) # there is a time effect 448
plmtest(ret1, c("time"), type=("bp")) # there is a time effect 449
plmtest(fet3, c("time"), type=("bp")) # there is a time effect 450
451 adf.test(panel3$recyc, k=2) 452
# Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 453
series is stationary 454
adf.test(panel3$fce, k=2) # stationary 455
adf.test(panel3$fce2, k=2) # stationary 456 adf.test(panel3$fce3, k=2) # stationary 457
458
# Model diagnostics 459
# All of these tests need to be performed individually for the models we would like to 460 compare
461 462
phtest(fet1,ret1) # fail to reject null (random) and state that the random effects are more 463
appropriate 464
phtest(fet3,ret3) # fixed effects are more appropriate 465 466
pcdtest(fet1) # not cross-sectional dependent 467
pcdtest(ret1) # not cross-sectional dependent 468
pcdtest(fet3) # not cross-sectional dependent 469
pcdtest(ret3) # not cross-sectional dependent 470
471 pbgtest(fet1) # there is serial correlation 472
pbgtest(fet3) # there is serial correlation 473
pbgtest(ret1) # there is serial correlation 474
pbgtest(ret3) # there is serial correlation 475 476
bptest(recyc ~ fce + factor(year), data = panel3, studentize=F) # there is heteroskedaticity 477
bptest(recyc ~ fce + fce2 + factor(year), data = panel3, studentize=F) # reject HO, thus 478
there is heteroskedaticity
479 bptest(recyc ~ fce + fce2 + fce3 + factor(year), data = panel3, studentize=F) # reject HO, 480
thus there is heteroskedaticity 481
482
#### Graphing against Regression Results ####
483
484 # Total 485
486
fun.1 <- function(x) 20.4389297 -3.4342419*x + 0.2043423*x^2 487
ggplot(data = panel1, aes(x=fce,y=total)) + 488
geom_point(aes(color=country)) + 489
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 490 Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 491
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 492
geom_function(aes(colour = "Cubic"), fun = fun.1) 493
494 # Landfill 495
496
fun.1 <- function(x) -115.942251 +27.431529*x - 1.543502*x^2 497
ggplot(data = panel2, aes(x=fce,y=land)) + 498
geom_point(aes(color=country)) +
499 labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 500
Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 501
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 502
geom_function(aes(colour = "Cubic"), fun = fun.1) 503
504 fun.1 <- function(x) 1.2224e+03 -4.1782e+02*x + 4.7747e+01*x^2 - 1.8161e+00*x^3 505
ggplot(data = panel2, aes(x=fce,y=land)) + 506
geom_point(aes(color=country)) + 507
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 508
Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 509 Calculation") + theme(legend.position = "", text = element_text(size = 20)) +
510
geom_function(aes(colour = "Cubic"), fun = fun.1) 511
512
# Recycling of Materials 513 514
fun.1 <- function(x) -1.0176e+03*x + 1.1106e+02*x^2 - -4.0295e+00*x^3 515
ggplot(data = panel3, aes(x=fce,y=recyc)) + 516
geom_point(aes(color=country)) + 517
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 518 Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 519
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 520
geom_function(aes(colour = "Cubic"), fun = fun.1) 521
522
fun.1 <- function(x) -13.606188 + 1.866839*x 523
ggplot(data = panel3, aes(x=fce,y=recyc)) + 524 geom_point(aes(color=country)) +
525
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 526
Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 527
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 528 geom_function(aes(colour = "Cubic"), fun = fun.1)
529 530
# Thesis Chapter Two 531
# Packaging Waste 532 533
require(readxl) 534
require(tidyverse) 535
require(ggpubr) 536
require(viridis) 537 require(plm) 538
require(lmtest) 539
require(tseries) 540
541
#### Descriptive Statistics ####
542
543 # All packaging waste generated for selected countries across years 544
packaging_waste_generated <- read_excel("Desktop/Indicators 545
Project/Descriptive/Transformed Data/packaging waste generated.xlsx") 546
547 all = packaging_waste_generated[,1:3] %>% rename(Region = Country) %>%
548
pivot_longer(-c(Region,Year), names_to = "type", values_to = "waste") 549
ggplot(data = all, aes(x=Year,y=waste,group=Region)) + 550
geom_point(aes(color=Region),size=4) + 551
geom_line(aes(color=Region)) +
552 labs(x="Year",y="Packaging Waste in kg per capita", title = "All Packaging Waste 553
Generated over Time", caption="Source: Eurostat") + 554
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 555
556
pwg = packaging_waste_generated %>% select(-"All packaging") %>% rename(Region = 557 Country) %>% pivot_longer(-c(Region,Year), names_to = "Type", values_to = "waste") 558
ggplot(data = pwg, aes(x=Year,y=waste,group=Region)) + 559
geom_point(aes(color=Region,shape=Type),size=3) + 560
geom_line(aes(color=Region)) + 561
facet_grid(.~Type,scales = "free") +
562 labs(x="Year",y="Packaging Waste in kg per capita", title = "Packaging Waste Generated 563
by Type over Time", caption="Source: Eurostat") + 564
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 565
566 # Total Recovery 567
total_recovery <- read_excel("Desktop/Indicators Project/Descriptive/Transformed 568
Data/total recovery.xlsx") 569
570
allry = total_recovery[,1:3] %>% rename(Region = Country)%>% pivot_longer(- 571 c(Region,Year), names_to = "type", values_to = "waste")
572
allryplot=ggplot(data = allry, aes(x=Year,y=waste,group=Region)) + 573
geom_point(aes(color=Region),size=4) + 574
geom_line(aes(color=Region)) + 575
labs(x="Year",y="Recovered Packaging Waste in kg per capita", title = "14A: All 576
Packaging Waste Recovered over Time", caption="Source: Eurostat") + 577 theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 578
579
# Percentage Version 580
pr <- read_excel("Desktop/Indicators Project/Descriptive/Transformed Data/percentage 581 recovered.xlsx")
582
prall = pr %>% filter(Type=="All") %>% pivot_longer(-c(Type,year), names_to = "Region", 583
values_to = "waste") 584
585 prallplot = ggplot(data = prall, aes(x=year,y=waste,group=Region)) + 586
geom_point(aes(color=Region),size=3) + 587
geom_line(aes(color=Region)) + 588
labs(x="Year",y="Percentage of Recovered Waste", title = "14B: All Recovered Packaging 589
Waste over Time", caption="Source: Eurostat") +
590 theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 591
592
ggarrange(allryplot,prallplot) 593
594
pwgry = total_recovery %>% select(-"All packaging") %>% rename(Region = Country) 595
%>% pivot_longer(-c(Region,Year), names_to = "Type", values_to = "waste") 596 ggplot(data = pwgry, aes(x=Year,y=waste,group=Region)) +
597
geom_point(aes(color=Region,shape=Type),size=3) + 598
geom_line(aes(color=Region)) + 599
facet_grid(.~Type,scales = "free") +
600 labs(x="Year",y="Recovered Packaging Waste in kg per capita", title = "Recovered 601
Packaging Waste by Type over Time", caption="Source: Eurostat") + 602
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 603
604
# PERCENTAGE VERSION
605 pro = pr %>% filter(Type!="All") %>% pivot_longer(-c(Type,year), names_to = "Region", 606
values_to = "waste") 607
608
ggplot(data = pro, aes(x=year,y=waste,group=Region)) + 609
geom_point(aes(color=Region,shape=Type),size=3) + 610 geom_line(aes(color=Region,linetype=Type)) + 611
facet_wrap(~Type,scales = "free_y",ncol=2) + 612
labs(x="Year",y="Percentage of Recovered Waste", title = "Recovered Packaging Waste 613
by Type over Time", caption="Source: Eurostat") + 614
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 615 616
# Total Recycling 617
618
total_recycling <- read_excel("Desktop/Indicators Project/Descriptive/Transformed 619 Data/total recycling.xlsx")
620 621
allre = total_recycling[,1:3] %>% rename(Region = Country)%>% pivot_longer(- 622
c(Region,Year), names_to = "type", values_to = "waste") 623
a = ggplot(data = allre, aes(x=Year,y=waste,group=Region)) + 624 geom_point(aes(color=Region),size=4) +
625
geom_line(aes(color=Region)) + 626
labs(x="Year",y="Recycled Packaging Waste in kg per capita", title = "16A: All Packaging 627
Waste Recycled over Time", caption="Source: Eurostat") + 628
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 629
630
pre <- read_excel("Desktop/Indicators Project/Descriptive/Transformed Data/percentage 631
recycled.xlsx") 632
preall = pre %>% filter(Type=="All") %>% pivot_longer(-c(Type,year), names_to = 633 "Region", values_to = "waste")
634 635
b = ggplot(data = preall, aes(x=year,y=waste,group=Region)) + 636
geom_point(aes(color=Region),size=4) + 637 geom_line(aes(color=Region)) +
638
labs(x="Year",y="Percentage of Recycled Waste", title = "16B: Recycled Packaging Waste 639
by Type over Time", caption="Source: Eurostat") + 640
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 641
642 ggarrange(a,b) 643
644
pwgre = total_recycling %>% select(-"All packaging") %>% rename(Region = Country) 645
%>% pivot_longer(-c(Region,Year), names_to = "Type", values_to = "waste") 646
ggplot(data = pwgre, aes(x=Year,y=waste,group=Region)) + 647
geom_point(aes(color=Region,shape=Type),size=3) + 648 geom_line(aes(color=Region)) +
649
facet_grid(.~Type,scales = "free") + 650
labs(x="Year",y="Recycled Packaging Waste in kg per capita", title = "Recycled Packaging 651
Waste by Type over Time", caption="Source: Eurostat") +
652 theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 653
654
# Percentage Version 655
656
preo = pre %>% filter(Type!="All") %>% pivot_longer(-c(Type,year), names_to = "Region", 657 values_to = "waste")
658 659
ggplot(data = preo, aes(x=year,y=waste,group=Region)) + 660
geom_point(aes(color=Region,shape=Type),size=3) + 661
geom_line(aes(color=Region,linetype=Type)) + 662 facet_wrap(~Type,scales = "free_y",ncol=2) + 663
labs(x="Year",y="Percentage of Recycled Waste", title = "Recycled Packaging Waste by 664
Type over Time", caption="Source: Eurostat") + 665
theme(legend.position = "bottom",axis.text.x = element_text(angle = 90)) 666
667 #### Preparing Data for Panel Regression ####
668 669
income1 <- read_excel("Desktop/Indicators Project/income.xlsx") 670
regions <- read_excel("Desktop/Indicators Project/regions.xlsx") 671 income = left_join(regions, income1, by = c("country"))
672 673
total <- read_excel("Desktop/Indicators Project/total.xlsx") 674
incomelong = income %>% pivot_longer(-c(country,region), names_to = "year", values_to = 675
"gdpo")
676 totallong = total %>% pivot_longer(-c(country), names_to = "year", values_to = "totalo") 677
totalw = left_join(totallong, incomelong, by = c("country","year")) 678
saveRDS(totalw, file = "totalw.rds") 679
680
paper <- read_excel("Desktop/Indicators Project/paper.xlsx") 681
paperlong = paper %>% pivot_longer(-c(country), names_to = "year", values_to = 682 "papero")
683
paperw = left_join(paperlong, incomelong, by = c("country","year")) 684
saveRDS(paperw, file = "paperw.rds") 685
686 plastic <- read_excel("Desktop/Indicators Project/plastic.xlsx") 687
plasticlong = plastic %>% pivot_longer(-c(country), names_to = "year", values_to = 688
"plastico") 689
plasticw = left_join(plasticlong, incomelong, by = c("country","year")) 690 saveRDS(plasticw, file = "plasticw.rds")
691 692
rm(list = ls()) 693
694
#### Panel Regression Descriptives ####
695 # plm::is.pbalanced(panel1$country,panel1$year) # check if it is a balanced panel 696
697
plm::is.pbalanced(panel1$country,panel1$year) # check if it is a balanced panel 698
totalw <- readRDS("~/Desktop/Indicators Project/totalw.rds") 699
summary(totalw) 700
701 panel1 = totalw %>% mutate(gdp = log(gdpo),gdp2 = 702
log(gdpo)*log(gdpo),gdp3=log(gdpo)*log(gdpo)*log(gdpo),total=log(totalo),.keep="unused"
703 ) 704
705 ggplot(data = panel1, aes(x=gdp,y=total)) + 706
geom_point(aes(color=country)) + 707
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Total Packaging 708
Waste Generated) in kg per capita", title = "Total Packaging Waste vs Income", 709
caption="Source: Eurostat") + theme(legend.position = "bottom") + 710 scale_color_viridis(discrete=TRUE)
711 712
ggplot(data = panel1, aes(x=gdp,y=total)) + 713
geom_point(aes(color=country)) + 714
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Total Packaging 715 Waste Generated) in kg per capita", title = "Total Packaging Waste vs Income by Country", 716
caption="Source: Eurostat") + theme(legend.position = "bottom") + 717
scale_color_viridis(discrete=TRUE) + 718
facet_wrap(~country,scales = "free") + 719
theme(legend.position = "none") 720 721
# PAPER 722
paperw <- readRDS("~/Desktop/Indicators Project/paperw.rds") 723
summary(paperw) 724 725
panel2 = paperw %>% mutate(gdp = log(gdpo),gdp2 = 726
log(gdpo)*log(gdpo),gdp3=log(gdpo)*log(gdpo)*log(gdpo),paper=log(papero),.keep="unus 727
ed") 728
729 ggplot(data = panel2, aes(x=gdp,y=paper)) + 730
geom_point(aes(color=country)) + 731
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Paper and 732
Cardboard Packaging Waste Generated) in kg per capita", title = "Paper and Cardboard 733
Packaging Waste vs Income", caption="Source: Eurostat") + theme(legend.position = 734
"bottom") + scale_color_viridis(discrete=TRUE) 735 736
ggplot(data = panel2, aes(x=gdp,y=paper)) + 737
geom_point(aes(color=country)) + 738
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Paper and 739 Cardboard Packaging Waste Generated) in kg per capita", title = "Paper and Cardboard 740
Packaging Waste vs Income by Country", caption="Source: Eurostat") + theme(legend.position 741
= "bottom") + scale_color_viridis(discrete=TRUE) + 742
facet_wrap(~country,scales = "free") + 743 theme(legend.position = "none")
744 745
# PLASTIC 746
plasticw <- readRDS("~/Desktop/Indicators Project/plasticw.rds") 747
summary(plasticw) 748 749
panel3 = plasticw %>% mutate(gdp = log(gdpo),gdp2 = 750
log(gdpo)*log(gdpo),gdp3=log(gdpo)*log(gdpo)*log(gdpo),plastic=log(plastico),.keep="unu 751
sed") 752
753
ggplot(data = panel3, aes(x=gdp,y=plastic)) + 754 geom_point(aes(color=country)) +
755
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 756
Waste Generated) in kg per capita", title = "Plastic Packaging Waste vs Income", 757
caption="Source: Eurostat") + theme(legend.position = "bottom") + 758 scale_color_viridis(discrete=TRUE)
759 760
ggplot(data = panel3, aes(x=gdp,y=plastic)) + 761
geom_point(aes(color=country)) + 762
labs(color="Country", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 763 Waste Generated) in kg per capita", title = "Plastic Packaging vs Income by Country", 764
caption="Source: Eurostat") + theme(legend.position = "bottom") + 765
scale_color_viridis(discrete=TRUE) + 766
facet_wrap(~country,scales = "free") + 767
theme(legend.position = "none") 768 769
#### Total Panel Regression PANEL1 ####
770 771
plm::is.pbalanced(panel1$country,panel1$year) # check if it is a balanced panel 772
773 # fixed effects 774
fe1 = plm(total ~ gdp + region, data = panel1,index=c("country","year"),model="within") # 775
total 776
summary(fe1, vcovHC) 777 778
# fixed effects with squared term 779
fe2 = plm(total ~ gdp + gdp2 + region, data = 780
panel1,index=c("country","year"),model="within") 781
summary(fe2, vcovHC) 782 783
# fixed effects with cubic term 784
fe3 = plm(total ~ gdp + gdp2 + gdp3 + region, data = 785
panel1,index=c("country","year"),model="within") 786
summary(fe3, vcovHC) 787
788 # random effects 789
re1 = plm(total ~ gdp + region, data = panel1,index=c("country","year"),model="random") 790
summary(re1, vcovHC) 791
792 #random effects with squared term 793
re2 = plm(total ~ gdp + gdp2 + region, data = 794
panel1,index=c("country","year"),model="random") 795
summary(re2, vcovHC) 796 797
#random effects with cubic term 798
re3 = plm(total ~ gdp + gdp2 + gdp3 + region, data = 799
panel1,index=c("country","year"),model="random") 800
summary(re3, vcovHC) 801 802
# fixed effects with time 803
fet1 = plm(total ~ gdp + factor(year) + region, data=panel1, model="within") 804
summary(fet1,vcovHC) 805
806
fet2 = plm(total ~ gdp + gdp2 + factor(year) + region, data=panel1, model="within") 807 summary(fet2,vcovHC) # with squared term.
808 809
fet3 = plm(total ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel1, 810
model="within")
811 summary(fet3,vcovHC) # with cubic term 812
813
# random effects with time 814
ret1 = plm(total ~ gdp + region + factor(year), data=panel1, 815
model="random",index=c("country","year")) 816 summary(ret1,vcovHC)
817 818
ret2 = plm(total ~ gdp + gdp2 + factor(year) + region, data=panel1, 819
model="random",index=c("country","year")) 820
summary(ret2,vcovHC) # with squared term 821 822
ret3 = plm(total ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel1, 823
model="random",index=c("country","year")) 824
summary(ret3,vcovHC) # with cubic term 825
826 # panel diagnostics 827
828
cor(panel1$gdp,panel1$total, method="pearson") # The correlation between the two 829
variables is 0.757032.
830 831
pt1 = plm(total ~ gdp + gdp2 + gdp3 + factor(year) + factor(country), data=panel1, 832
model="pooling") 833
plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 834
pooled OLS is not appropriate
835 pt2 = plm(total ~ gdp + gdp2 + factor(year) + factor(country), data=panel1, 836
model="pooling") 837
plmtest(pt2, type=c("bp")) # there is a panel effect 838
pt3 = plm(total ~ gdp + factor(year) + factor(country), data=panel1, model="pooling") 839
plmtest(pt3, type=c("bp")) # there is a panel effect 840
841 pFtest(fet1,fe1) # reject HO, and state there is a time effect 842
pFtest(ret1,re1) # there is a time effect 843
pFtest(ret3,re3) # there is a time effect 844
plmtest(fe1, c("time"), type=("bp")) # do not reject HO and state there is no time effect 845 plmtest(re1, c("time"), type=("bp")) # there is no time effect
846 847
adf.test(panel1$total, k=2) 848
# Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 849 series is stationary
850
adf.test(panel1$gdp, k=2) 851
852
# Model diagnostics 853
# All of these tests need to be performed individually for the models we would like to 854 compare
855 856
phtest(fet1,ret1) # fail to reject null (random) and state that the random effects is more 857
appropriate 858
phtest(fet3,ret3) # random effects is more appropriate 859
phtest(fe1,re1) # random effects is more appropriate 860 phtest(fe3,re3) # random effects is more appropriate 861
862
pcdtest(re1) # is cross-sectional dependent 863
pcdtest(re3) # is cross-sectional dependent 864 pcdtest(ret1) # not cross-sectional dependent 865
866
pbgtest(re1) # there is serial correlation 867
pbgtest(re3) # there is serial correlation 868
pbgtest(ret1) # there is serial correlation 869 870
bptest(total ~ gdp + factor(year), data = panel1, studentize=F) # fail to reject HO, thus 871
there's no heteroskedaticity 872
bptest(total ~ gdp + gdp2 + gdp3 + factor(year), data = panel1, studentize=F) # reject 873
HO, thus there is heteroskedaticity 874 875
876
#### Paper Panel Regression PANEL 2####
877 878
plm::is.pbalanced(panel2$country,panel2$year) # check if it is a balanced panel 879 880
# fixed effects 881
fe1 = plm(paper ~ gdp + region, data = panel2,index=c("country","year"),model="within") 882
# paper
883 summary(fe1, vcovHC) 884
885
# fixed effects with squared term 886
fe2 = plm(paper ~ gdp + gdp2 + region, data = 887
panel2,index=c("country","year"),model="within") 888 summary(fe2, vcovHC)
889 890
# fixed effects with cubic term 891
fe3 = plm(paper ~ gdp + gdp2 + gdp3 + region, data = 892
panel2,index=c("country","year"),model="within") 893
summary(fe3, vcovHC) 894 895
# random effects 896
re1 = plm(paper ~ gdp + region, data = 897
panel2,index=c("country","year"),model="random") 898 summary(re1, vcovHC)
899 900
#random effects with squared term 901
re2 = plm(paper ~ gdp + gdp2 + region, data = 902 panel2,index=c("country","year"),model="random") 903
summary(re2, vcovHC) 904
905
#random effects with cubic term 906
re3 = plm(paper ~ gdp + gdp2 + gdp3 + region, data = 907 panel2,index=c("country","year"),model="random")
908
summary(re3, vcovHC) 909
910
# fixed effects with time 911
fet1 = plm(paper ~ gdp + factor(year) + region, data=panel2, model="within") 912
summary(fet1,vcovHC) 913 914
fet2 = plm(paper ~ gdp + gdp2 + factor(year) + region, data=panel2, model="within") 915
summary(fet2,vcovHC) # with squared term.
916
917 fet3 = plm(paper ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel2, 918
model="within") 919
summary(fet3,vcovHC) # with cubic term 920
921
# random effects with time
922 ret1 = plm(paper ~ gdp + region + factor(year), data=panel2, 923
model="random",index=c("country","year")) 924
summary(ret1,vcovHC) 925
926
ret2 = plm(paper ~ gdp + gdp2 + factor(year) + region, data=panel2, 927 model="random",index=c("country","year"))
928
summary(ret2,vcovHC) # with squared term 929
930
ret3 = plm(paper ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel2, 931
model="random",index=c("country","year")) 932 summary(ret3,vcovHC) # with cubic term 933
934
# panel diagnostics 935
936 cor(panel2$gdp,panel2$paper, method="pearson") # The correlation between the two 937
variables is 0.757032.
938 939
pt1 = plm(paper ~ gdp + gdp2 + gdp3 + factor(year) + factor(country), data=panel2, 940
model="pooling")
941 plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 942
pooled OLS is not appropriate 943
pt2 = plm(paper ~ gdp + gdp2 + factor(year) + factor(country), data=panel2, 944
model="pooling") 945
plmtest(pt2, type=c("bp")) # there is a panel effect 946
pt3 = plm(paper ~ gdp + factor(year) + factor(country), data=panel2, model="pooling") 947 plmtest(pt3, type=c("bp")) # there is a panel effect
948
949
pFtest(fet1,fe1) # reject HO, and state there is a time effect 950
pFtest(ret1,re1) # there is a time effect
951 plmtest(fe1, c("time"), type=("bp")) # do not reject HO and state there is no time effect 952
plmtest(re1, c("time"), type=("bp")) # there is no time effect 953
954
adf.test(panel2$paper, k=2)
955 # Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 956
series is stationary 957
adf.test(panel2$gdp, k=2) 958
959
# Model diagnostics
960 # All of these tests need to be performed individually for the models we would like to 961
compare 962
963
phtest(fet1,ret1) # fail to reject null (random) and state that the random effects is more 964
appropriate 965
phtest(fe1,re1) # random effects is more appropriate 966 967
pcdtest(re1) # is cross-sectional dependent 968
pcdtest(ret1) # is not cross-sectional dependent 969
970 pbgtest(re1) # there is serial correlation 971
pbgtest(ret1) # there is serial correlation 972
973
bptest(paper ~ gdp + factor(year), data = panel2, studentize=F) # fail to reject HO, thus 974
there's no heteroskedaticity 975 976
#### Plastic Panel Regression PANEL 3####
977 978
plm::is.pbalanced(panel3$country,panel3$year) # check if it is a balanced panel 979
980 # fixed effects 981
fe1 = plm(plastic ~ gdp + region, data = panel3,index=c("country","year"),model="within") 982
# plastic 983
summary(fe1, vcovHC) 984
985 # fixed effects with squared term 986
fe2 = plm(plastic ~ gdp + gdp2 + region, data = 987
panel3,index=c("country","year"),model="within") 988
summary(fe2, vcovHC) 989 990
# fixed effects with cubic term 991
fe3 = plm(plastic ~ gdp + gdp2 + gdp3 + region, data = 992
panel3,index=c("country","year"),model="within") 993
summary(fe3, vcovHC) 994 995
# random effects 996
re1 = plm(plastic ~ gdp + region, data = 997
panel3,index=c("country","year"),model="random") 998
summary(re1, vcovHC) 999
1000
#random effects with squared term 1001
re2 = plm(plastic ~ gdp + gdp2 + region, data = 1002
panel3,index=c("country","year"),model="random") 1003
summary(re2, vcovHC) 1004
1005
#random effects with cubic term 1006
re3 = plm(plastic ~ gdp + gdp2 + gdp3 + region, data = 1007
panel3,index=c("country","year"),model="random") 1008
summary(re3, vcovHC) 1009
1010
# fixed effects with time 1011
fet1 = plm(plastic ~ gdp + factor(year) + region, data=panel3, model="within") 1012
summary(fet1,vcovHC) 1013
1014
fet2 = plm(plastic ~ gdp + gdp2 + factor(year) + region, data=panel3, model="within") 1015
summary(fet2,vcovHC) # with squared term.
1016 1017
fet3 = plm(plastic ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel3, 1018
model="within") 1019
summary(fet3,vcovHC) # with cubic term 1020
1021
# random effects with time 1022
ret1 = plm(plastic ~ gdp + region + factor(year), data=panel3, 1023
model="random",index=c("country","year")) 1024
summary(ret1,vcovHC) 1025
1026
ret2 = plm(plastic ~ gdp + gdp2 + factor(year) + region, data=panel3, 1027
model="random",index=c("country","year")) 1028
summary(ret2,vcovHC) # with squared term 1029
1030
ret3 = plm(plastic ~ gdp + gdp2 + gdp3 + factor(year) + region, data=panel3, 1031
model="random",index=c("country","year")) 1032
summary(ret3,vcovHC) # with cubic term 1033
1034
# panel diagnostics 1035
1036
cor(panel3$gdp,panel3$plastic, method="pearson") # The correlation between the two 1037
variables is 0.757032.
1038 1039
pt1 = plm(plastic ~ gdp + gdp2 + gdp3 + factor(year) + factor(country), data=panel3, 1040
model="pooling") 1041
plmtest(pt1, type=c("bp")) # reject HO, and state that there is a panel effect, therefore 1042
pooled OLS is not appropriate 1043
pt2 = plm(plastic ~ gdp + gdp2 + factor(year) + factor(country), data=panel3, 1044
model="pooling") 1045
plmtest(pt2, type=c("bp")) # there is a panel effect 1046
pt3 = plm(plastic ~ gdp + factor(year) + factor(country), data=panel3, model="pooling") 1047
plmtest(pt3, type=c("bp")) # there is a panel effect 1048
1049
pFtest(fet1,fe1) # no time effect 1050
pFtest(ret1,re1) # no time effect 1051
pFtest(fet3,fe3) # no time effect 1052
pFtest(ret3,re3) # no time effect 1053
plmtest(fe3, c("time"), type=("bp")) # do not reject HO and state there is no time effect 1054
plmtest(re3, c("time"), type=("bp")) # there is no time effect 1055
1056
adf.test(panel3$plastic, k=2) 1057
# Because p-value < 0.05, therefore the series dos not have unit roots. This means that the 1058
series is stationary 1059
adf.test(panel3$gdp, k=2) 1060
1061
# Model diagnostics 1062
# All of these tests need to be performed individually for the models we would like to 1063
compare 1064
1065
phtest(fe1,re1) # random effects is more appropriate 1066
phtest(fe3,re3) # fixed effects is more appropriate 1067
1068
pcdtest(re1) # is not cross-sectional dependent 1069
pcdtest(fe3) # is not cross-sectional dependent 1070
1071
pbgtest(re1) # there is serial correlation 1072
pbgtest(fe3) # there is serial correlation 1073
1074
bptest(plastic ~ gdp + gdp2 + gdp3 + factor(year), data = panel3, studentize=F) # fail to 1075
reject HO, thus there's no heteroskedaticity 1076
bptest(plastic ~ gdp + factor(year), data = panel3, studentize=F) # no heteroskedaticity 1077
1078
# Graphing model 27 and 28 1079
fun.1 <- function(x) -101.018596*x + 10.131925*x^2 -0.335831*x^3 1080
cubic = ggplot(data = panel3, aes(x=gdp,y=plastic)) + 1081
geom_point(aes(color=country)) + 1082
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 1083
Waste Generated) in kg per capita", title = "Cubic Result", caption="Source: Eurostat/Own 1084
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 1085
geom_function(aes(colour = "Cubic"), fun = fun.1) 1086
1087
fun.2 <- function(x) -3.331390 + 0.643917*x 1088
linear = ggplot(data = panel3, aes(x=gdp,y=plastic)) + 1089
geom_point(aes(color=country)) + 1090
labs(color="Key", x="ln(GDP) in current prices PPS per capita",y="ln(Plastic Packaging 1091
Waste Generated) in kg per capita", title = "Linear Result", caption="Source: Eurostat/Own 1092
Calculation") + theme(legend.position = "", text = element_text(size = 20)) + 1093
geom_function(aes(colour = "Linear"), fun = fun.2) 1094
1095
ggarrange(cubic,linear) 1096
1097
# Comparing correlation 1098
cor(paperw$papero,paperw$gdpo, method="pearson") #0.6305646 1099
cor(panel2$gdp,panel2$paper,method="pearson") #0.7649613 1100
1101
cor(plasticw$plastico,plasticw$gdpo, method="pearson") #0.6538363 1102
cor(panel3$gdp,panel3$plastic,method="pearson") #0.7090483 1103
1104
cor(totalw$totalo,totalw$gdpo,method="pearson") #0.6890545 1105
cor(panel1$gdp,panel1$total,method="pearson") #0.757032 1106
1107
# Thesis Chapter Three 1108
# Waste Trade Network Analysis 1109
1110
# Libraries 1111
library(readr) 1112
library(tidyverse) 1113
library(igraph) 1114
library(readxl) 1115
library(countrycode) 1116
library(RColorBrewer) 1117
library(plotrix) 1118
library(viridis) 1119
library(ggpubr) 1120
library(ggsci) 1121
library(cowplot) 1122
set.seed(1234) 1123
1124
#### Descriptive Statistics ####
1125 1126
summary(import_2020$QUANTITY_IN_100KG) 1127
summary(export_2020$QUANTITY_IN_100KG) 1128
1129
euim = import_2020 %>% group_by(REPORTER) %>% summarise(import = 1130
sum(QUANTITY_IN_100KG)) 1131
euex = export_2020 %>% group_by(REPORTER) %>% summarise(export = 1132
sum(QUANTITY_IN_100KG)) 1133
euexpart = export_2020 %>% group_by(PARTNER) %>% summarise(export = 1134
sum(QUANTITY_IN_100KG)) 1135
df = euim %>% left_join(euex,by="REPORTER") %>%
1136
mutate(reporter = word(REPORTER,1)) %>%
1137
# mutate(net.export = export - import) %>%
1138
select(-c(REPORTER)) %>%
1139
rename(Import = import,Export = export) %>% pivot_longer(-c(reporter), names_to = 1140
"Type", values_to = "quantity") %>%
1141
mutate(tons = quantity * 100 / 1000) 1142
1143
ggplot(df, aes(fill=Type, y=tons, x=reporter)) + 1144
geom_bar(position="dodge", stat="identity") + 1145
labs(color="", x="Country",y="Paper and Cardboard Waste in metric tons", title = "Trade 1146
of Waste Paper and Cardboard by Country and Type ", caption="Source: Eurostat") + 1147
theme(legend.position = "bottom",axis.text.x=element_text(angle = 90)) 1148
1149
#### Data Preparation ####
1150 1151
# IMPORT NETWORK 1152
import_2020 <- read_excel("Downloads/import 2020.xlsx") 1153
import = import_2020 %>% mutate(reporter = word(REPORTER,1), partner = 1154
word(PARTNER,1)) %>%
1155
mutate_at("partner", str_replace, "Viet", "Vietnam") %>%
1156
rename(quantity = QUANTITY_IN_100KG) %>%
1157
select(c(reporter,partner,quantity)) %>%
1158
filter(quantity >= median(quantity),) %>%
1159
mutate_at("partner", str_replace, "United", "United States") %>%
1160
