workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                      ##
########################################################################################
library(xtable)
library(gdata)
library(dplyr)
library(reshape2)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(bindrcpp)
########################################################################################
## Load datasets                                                                      ##
########################################################################################
## GTAP bilateral margins and trade  (V 9.0)
vtwr.9.0 <- read.csv("../DATA/RawData/other/VTWR_9aY11.CSV")## sol( solution.dir = "../DATA/GTAP_Data/", solution.out = "VTWR_9aY11.HAR", csv.out = "VTWR_9aY11.CSV", header = "VTWR")
vxmd.9.0 <- read.csv("../DATA/RawData/other/VXMD_9aY11.CSV")
## Our new estimated shares:
nuvi <- read.csv("../DATA/output/OURSHARES.csv")
gtap.names <- read.csv("../DATA/RawData/other/gtapCodesComplete.csv")
bildist <- read.csv("../DATA/RawData/other/dist_cepii.csv") %>%
select( iso_o, iso_d, dist) %>%
dplyr::rename( o = iso_o, d = iso_d ) %>%
mutate( o = tolower(o),  d = tolower(d))
########################################################################################
## Clean and merge datasets                                                           ##
########################################################################################
vtwr <- rename.vars( data = vtwr.9.0,
from = c("MARG_COMM", "TRAD_COMM", "REG", "REG.1", "Value"),
to = c( "m", "g", "o", "d", "vtwr") )
vxmd <- rename.vars( data = vxmd.9.0,
from = c("TRAD_COMM", "REG", "REG.1", "Value"),
to = c("g", "o", "d", "vxmd") )
nuvi.mod <- nuvi %>%
mutate( o = tolower(o),
d = tolower(d),
m = recode( m , "Water" = "wtp", "Other" = "otp", "Air" =  "atp"),
value_MillUS = value_usd/1e6,
weight_ton = weight_kg/1e6) %>%
left_join( . , gtap.names[, c("g", "Code")],
by = ( "g") ) %>%
select( o, d, Code, m, value_MillUS, weight_ton, share, wv_ratio ) %>%
dplyr::rename( g = Code, wv_ratio.nuvi = wv_ratio, modal.share.nuvi = share )
## GTAP sectors to exclude: services, but also sugar cane and raw
## milk, where trade is marginal:
exclude <- c("c_b", "rmk", "ely","gdt","wtr","cns","trd","otp","wtp","atp","cmn","ofi","isr","obs","ros","osg","dwe")
## Weight-to-value ratios:
wv.ratios.by.g <- nuvi.mod %>%
group_by(o, d, g) %>%
summarise(wv.ratio = mean(wv_ratio.nuvi))
## Create a dataset that is square, that is, it has all the possible
## combinations of exporters, importers, commodities, and modes.
## Joining the other datasets to this dataset should help to track
## missing values while ensuring that all possible modal shares,
## including when they are zero, are used in the analysis:
origins <- unique(vxmd$o) ## 141
destinations <- unique(vxmd$d) ## 141
sectors <- setdiff( unique(vxmd$g), exclude) ## 42 sectors
modes <- unique(vtwr$m) ## 3 modes
grid <- expand.grid( origins, destinations, sectors, modes) %>%
dplyr::rename( o = Var1, d = Var2, g = Var3, m = Var4 ) %>%
filter( o != d ) %>%
arrange(o , d, g, m)
## Add trade (the NAS here is because there is no trade: set them to
## zero):
dat <- left_join(grid, vxmd, by = c("o", "d", "g") ) %>%
left_join(. , vtwr, by = c("o", "d", "g", "m") )
dat[is.na(dat)] <- 0
## To calculate modal shares we need to add margins across modes,
## these are then merged back into dat, and the ratios (or modal
## shares) are taken:
transport.margins <- dat %>%
group_by(o, d, g)  %>%
summarise(transport.margin.gtap = sum(vtwr, na.rm = TRUE))
## Calculate modal shares:
dat <- left_join(dat, transport.margins, by = c("o", "d", "g") ) %>%
mutate(modal.share.gtap = vtwr/transport.margin.gtap)
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                    ##
########################################################################################
library(plyr)
library(dplyr)
library(gdata)
library(reshape2)
library(tidyr)
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                    ##
########################################################################################
library(plyr)
library(dplyr)
library(gdata)
library(reshape2)
library(tidyr)
########################################################################################
## Load Database                                                                      ##
########################################################################################
# Load database with the information to be used.
load("../DATA/my_data/DF.RData")
DF2 <- DF
rm(DF)
# The equations to be estimated in do.r are:
# Sea transport: ln(f/v)_{ckt} = -4.16 + 0.393 ln(w/v) + 0.351 ln(oil) + 0.027 DIST + 0.0063 ln(oil)*DIST
# Air transport: ln(f/v)_{ckt} = -3.80 + 0.435 ln(w/v) + 0.209 ln(jetfuel) + 0.033 DIST + 0.018 ln(jetfuel)*DIST
# Road transport: ln(f/v)_{ckt} = -8.18 + 0.234 ln(w/v) + 0.862 ln(oil) + 1.21 ln(DIST) - 0.373ln(oil)*ln(DIST)
# Rail transport: ln(f/v)_{ckt} = -4.37 + 0.54 ln(w/v) + 0.079 ln(oil) - 0.90 ln(DIST) - 0.224 ln(oil)*ln(DIST)
# from: "The Wage Effects of Offshoring: Evidence from Danish Matched Worker-Firm Data" Hummels et al. 2013
DF2[,c('fv')] <- NA
DF2 <- within(DF2, fv[DF2$transp_mode == 'Water'] <- -4.16
+ 0.393*log(wv_ratio[DF2$transp_mode == 'Water'])
+ 0.351*log(oil[DF2$transp_mode == 'Water'])
+ 0.027*dist[DF2$transp_mode == 'Water']
+ 0.0063*log(oil[DF2$transp_mode == 'Water'])*(dist[DF2$transp_mode == 'Water']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Air'] <-  -3.80
+ 0.435*log(wv_ratio[DF2$transp_mode == 'Air'])
+ 0.209*log(jet[DF2$transp_mode == 'Air'])
+ 0.033*log(dist[DF2$transp_mode == 'Air'])
+ 0.018*log(jet[DF2$transp_mode == 'Air'])*(dist[DF2$transp_mode == 'Air']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Road'] <- -8.18
+ 0.234*log(wv_ratio[DF2$transp_mode == 'Road'])
+ 0.862*log(oil[DF2$transp_mode == 'Road'])
+ 1.21*log(dist[DF2$transp_mode == 'Road'])
- 0.373*log(oil[DF2$transp_mode == 'Road'])*log(dist[DF2$transp_mode == 'Road']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Rail'] <- -4.37
+ 0.54*log(wv_ratio[DF2$transp_mode == 'Rail'])
+ 0.079*log(oil[DF2$transp_mode == 'Rail'])
- 0.90*log(dist[DF2$transp_mode == 'Rail'])
- 0.224*log(oil[DF2$transp_mode == 'Rail'])*log(dist[DF2$transp_mode == 'Rail']))
# Getting Tau
DF2[,c('tau')] <- exp(DF2$fv)
#ad-valorem transport margins are CIF/FOB. If origin price is p and delivered price
#including margins is p*, and the advalorem margin is tau, then CIF/FOB gives you p*/p =
#(1+tau). CIF/FOB is in the column called "cif_fob"
# Getting the missing value, cif or fob. If the reported flow is cif, we need to find fob; if the reported flow is cif, we need
# to obtain cif. Therefore, if the reported values is cif, value2 will be the fob value and the other way around.
# cif = fob*(1+tau)
# fob = cif/(1+tau)
DF2[,c('value2')] <- NA
DF2 <- within(DF2, value2[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif']/(1 + tau[DF2$flow == 'cif']))
DF2 <- within(DF2, value2[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob']*(1 + tau[DF2$flow == 'fob']))
# We define a margin as cif - fob
DF2[,c('margin')] <- NA
DF2 <- within(DF2, margin[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'] - value2[DF2$flow == 'cif'])
DF2 <- within(DF2, margin[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'] - value_usd[DF2$flow == 'fob'])
# Getting the margin shares: margin/[cif_{total} - fob_{total}]
DF2[,c('cif_total')] <- 0
DF2[,c('fob_total')] <- 0
DF2 <- within(DF2, cif_total[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'cif'] <- value2[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob'])
DF2 <- within(DF2, cif_total[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'])
DF2[,c('odg')] <- paste(DF2$ISO_orig, DF2$ISO_dest, DF2$GTAP, sep="_")
DF2$cif_total <- ave(DF2$cif_total, DF2$odg, FUN=sum)
DF2$fob_total <- ave(DF2$fob_total, DF2$odg, FUN=sum)
DF2$share <- DF2$margin/(DF2$cif_total - DF2$fob_total)
#NANS <- DF2[is.na(DF2$share),]
DF2 <- DF2[!(is.na(DF2$share)),]
DF2 <- DF2[order(DF2$odg),]
row.names(DF2) <- NULL
# Column DF2$share is the one containing our margin shares. This one includes FOUR modes (air, sea, rail and road)
DT <- subset(DF2, select=c("ISO_orig", "ISO_dest", "value_usd", "weight_kg", "dist","GTAP", "transp_mode", "share", "wv_ratio"))
# We will repoort shares for air, water, and other (rail+ road)
DT[,c('mode')] <- NA
DT <- within(DT, mode[DT$transp_mode == "Rail"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Road"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Water"] <- 'Water')
DT <- within(DT, mode[DT$transp_mode == "Air"] <- 'Air')
DT$transp_mode <- NULL
DT[,c('odgm')] <- paste(DT$ISO_orig, DT$ISO_dest, DT$GTAP, DT$mode, sep="_")
DT$share <- ave(DT$share, DT$odgm, FUN=sum)
# "Other" is duplicated because it appears once where the mode was "rail" and once more where the mode was "road".
DT <- DT[!duplicated(DT$odgm), ]
DT <- DT[order(DT$odg),]
DT <- plyr::rename(DT, c("ISO_orig"="o", "ISO_dest"="d", "GTAP"="g", "mode"="m", "value_usd"="value_usd", "weight_kg"="weight_kg", "dist"="dist",
"share"="share", "wv_ratio"="wv_ratio", "odgm"="odgm" ))
write.csv(DT, "../DATA/output/OURSHARES.csv")
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
## Load R packages                                                                      ##
########################################################################################
library(xtable)
library(gdata)
library(dplyr)
library(reshape2)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(bindrcpp)
#################################
vtwr.9.0 <- read.csv("../DATA/RawData/other/VTWR_9aY11.CSV")## sol( solution.dir = "../DATA/GTAP_Data/", solution.out = "VTWR_9aY11.HAR", csv.out = "VTWR_9aY11.CSV", header = "VTWR")
vxmd.9.0 <- read.csv("../DATA/RawData/other/VXMD_9aY11.CSV")
nuvi <- read.csv("../DATA/output/OURSHARES.csv")
gtap.names <- read.csv("../DATA/RawData/other/gtapCodesComplete.csv")
bildist <- read.csv("../DATA/RawData/other/dist_cepii.csv") %>%
select( iso_o, iso_d, dist) %>%
dplyr::rename( o = iso_o, d = iso_d ) %>%
mutate( o = tolower(o),  d = tolower(d))
########################################################################################
## Clean and merge datasets                                                           ##
########################################################################################
vtwr <- rename.vars( data = vtwr.9.0,
from = c("MARG_COMM", "TRAD_COMM", "REG", "REG.1", "Value"),
to = c( "m", "g", "o", "d", "vtwr") )
vxmd <- rename.vars( data = vxmd.9.0,
from = c("TRAD_COMM", "REG", "REG.1", "Value"),
to = c("g", "o", "d", "vxmd") )
nuvi.mod <- nuvi %>%
mutate( o = tolower(o),
d = tolower(d),
m = recode( m , "Water" = "wtp", "Other" = "otp", "Air" =  "atp"),
value_MillUS = value_usd/1e6,
weight_ton = weight_kg/1e6) %>%
left_join( . , gtap.names[, c("g", "Code")],
by = ( "g") ) %>%
select( o, d, Code, m, value_MillUS, weight_ton, share, wv_ratio ) %>%
dplyr::rename( g = Code, wv_ratio.nuvi = wv_ratio, modal.share.nuvi = share )
## GTAP sectors to exclude: services, but also sugar cane and raw
## milk, where trade is marginal:
exclude <- c("c_b", "rmk", "ely","gdt","wtr","cns","trd","otp","wtp","atp","cmn","ofi","isr","obs","ros","osg","dwe")
## Weight-to-value ratios:
wv.ratios.by.g <- nuvi.mod %>%
group_by(o, d, g) %>%
summarise(wv.ratio = mean(wv_ratio.nuvi))
origins <- unique(vxmd$o) ## 141
destinations <- unique(vxmd$d) ## 141
sectors <- setdiff( unique(vxmd$g), exclude) ## 42 sectors
modes <- unique(vtwr$m) ## 3 modes
grid <- expand.grid( origins, destinations, sectors, modes) %>%
dplyr::rename( o = Var1, d = Var2, g = Var3, m = Var4 ) %>%
filter( o != d ) %>%
arrange(o , d, g, m)
## Add trade (the NAS here is because there is no trade: set them to
## zero):
dat <- left_join(grid, vxmd, by = c("o", "d", "g") ) %>%
left_join(. , vtwr, by = c("o", "d", "g", "m") )
dat[is.na(dat)] <- 0
## To calculate modal shares we need to add margins across modes,
## these are then merged back into dat, and the ratios (or modal
## shares) are taken:
transport.margins <- dat %>%
group_by(o, d, g)  %>%
summarise(transport.margin.gtap = sum(vtwr, na.rm = TRUE))
## Calculate modal shares:
dat <- left_join(dat, transport.margins, by = c("o", "d", "g") ) %>%
mutate(modal.share.gtap = vtwr/transport.margin.gtap)
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                    ##
########################################################################################
library(plyr)
library(dplyr)
library(gdata)
library(reshape2)
library(tidyr)
########################################################################################
## Load Database                                                                      ##
########################################################################################
# Load database with the information to be used.
load("../DATA/my_data/DF.RData")
DF2 <- DF
rm(DF)
# The equations to be estimated in do.r are:
# Sea transport: ln(f/v)_{ckt} = -4.16 + 0.393 ln(w/v) + 0.351 ln(oil) + 0.027 DIST + 0.0063 ln(oil)*DIST
# Air transport: ln(f/v)_{ckt} = -3.80 + 0.435 ln(w/v) + 0.209 ln(jetfuel) + 0.033 DIST + 0.018 ln(jetfuel)*DIST
# Road transport: ln(f/v)_{ckt} = -8.18 + 0.234 ln(w/v) + 0.862 ln(oil) + 1.21 ln(DIST) - 0.373ln(oil)*ln(DIST)
# Rail transport: ln(f/v)_{ckt} = -4.37 + 0.54 ln(w/v) + 0.079 ln(oil) - 0.90 ln(DIST) - 0.224 ln(oil)*ln(DIST)
# from: "The Wage Effects of Offshoring: Evidence from Danish Matched Worker-Firm Data" Hummels et al. 2013
DF2[,c('fv')] <- NA
DF2 <- within(DF2, fv[DF2$transp_mode == 'Water'] <- -4.16
+ 0.393*log(wv_ratio[DF2$transp_mode == 'Water'])
+ 0.351*log(oil[DF2$transp_mode == 'Water'])
+ 0.027*dist[DF2$transp_mode == 'Water']
+ 0.0063*log(oil[DF2$transp_mode == 'Water'])*(dist[DF2$transp_mode == 'Water']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Air'] <-  -3.80
+ 0.435*log(wv_ratio[DF2$transp_mode == 'Air'])
+ 0.209*log(jet[DF2$transp_mode == 'Air'])
+ 0.033*log(dist[DF2$transp_mode == 'Air'])
+ 0.018*log(jet[DF2$transp_mode == 'Air'])*(dist[DF2$transp_mode == 'Air']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Road'] <- -8.18
+ 0.234*log(wv_ratio[DF2$transp_mode == 'Road'])
+ 0.862*log(oil[DF2$transp_mode == 'Road'])
+ 1.21*log(dist[DF2$transp_mode == 'Road'])
- 0.373*log(oil[DF2$transp_mode == 'Road'])*log(dist[DF2$transp_mode == 'Road']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Rail'] <- -4.37
+ 0.54*log(wv_ratio[DF2$transp_mode == 'Rail'])
+ 0.079*log(oil[DF2$transp_mode == 'Rail'])
- 0.90*log(dist[DF2$transp_mode == 'Rail'])
- 0.224*log(oil[DF2$transp_mode == 'Rail'])*log(dist[DF2$transp_mode == 'Rail']))
# Getting Tau
DF2[,c('tau')] <- exp(DF2$fv)
#ad-valorem transport margins are CIF/FOB. If origin price is p and delivered price
#including margins is p*, and the advalorem margin is tau, then CIF/FOB gives you p*/p =
#(1+tau). CIF/FOB is in the column called "cif_fob"
# Getting the missing value, cif or fob. If the reported flow is cif, we need to find fob; if the reported flow is cif, we need
# to obtain cif. Therefore, if the reported values is cif, value2 will be the fob value and the other way around.
# cif = fob*(1+tau)
# fob = cif/(1+tau)
DF2[,c('value2')] <- NA
DF2 <- within(DF2, value2[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif']/(1 + tau[DF2$flow == 'cif']))
DF2 <- within(DF2, value2[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob']*(1 + tau[DF2$flow == 'fob']))
# We define a margin as cif - fob
DF2[,c('margin')] <- NA
DF2 <- within(DF2, margin[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'] - value2[DF2$flow == 'cif'])
DF2 <- within(DF2, margin[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'] - value_usd[DF2$flow == 'fob'])
# Getting the margin shares: margin/[cif_{total} - fob_{total}]
DF2[,c('cif_total')] <- 0
DF2[,c('fob_total')] <- 0
DF2 <- within(DF2, cif_total[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'cif'] <- value2[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob'])
DF2 <- within(DF2, cif_total[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'])
DF2[,c('odg')] <- paste(DF2$ISO_orig, DF2$ISO_dest, DF2$GTAP, sep="_")
DF2$cif_total <- ave(DF2$cif_total, DF2$odg, FUN=sum)
DF2$fob_total <- ave(DF2$fob_total, DF2$odg, FUN=sum)
DF2$share <- DF2$margin/(DF2$cif_total - DF2$fob_total)
#NANS <- DF2[is.na(DF2$share),]
DF2 <- DF2[!(is.na(DF2$share)),]
DF2 <- DF2[order(DF2$odg),]
row.names(DF2) <- NULL
# Column DF2$share is the one containing our margin shares. This one includes FOUR modes (air, sea, rail and road)
DT <- subset(DF2, select=c("ISO_orig", "ISO_dest", "value_usd", "weight_kg", "dist","GTAP", "transp_mode", "share", "wv_ratio"))
# We will repoort shares for air, water, and other (rail+ road)
DT[,c('mode')] <- NA
DT <- within(DT, mode[DT$transp_mode == "Rail"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Road"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Water"] <- 'Water')
DT <- within(DT, mode[DT$transp_mode == "Air"] <- 'Air')
DT$transp_mode <- NULL
DT[,c('odgm')] <- paste(DT$ISO_orig, DT$ISO_dest, DT$GTAP, DT$mode, sep="_")
DT$share <- ave(DT$share, DT$odgm, FUN=sum)
# "Other" is duplicated because it appears once where the mode was "rail" and once more where the mode was "road".
DT <- DT[!duplicated(DT$odgm), ]
DT <- DT[order(DT$odg),]
DT <- plyr::rename(DT, c("ISO_orig"="o", "ISO_dest"="d", "GTAP"="g", "mode"="m", "value_usd"="value_usd", "weight_kg"="weight_kg", "dist"="dist",
"share"="share", "wv_ratio"="wv_ratio", "odgm"="odgm" ))
write.csv(DT, "../DATA/output/OURSHARES.csv")
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                    ##
########################################################################################
library(plyr)
library(dplyr)
library(gdata)
library(reshape2)
library(tidyr)
########################################################################################
## Load Database                                                                      ##
########################################################################################
# Load database with the information to be used.
load("../DATA/my_data/DF.RData")
DF2 <- DF
rm(DF)
# The equations to be estimated in do.r are:
# Sea transport: ln(f/v)_{ckt} = -4.16 + 0.393 ln(w/v) + 0.351 ln(oil) + 0.027 DIST + 0.0063 ln(oil)*DIST
# Air transport: ln(f/v)_{ckt} = -3.80 + 0.435 ln(w/v) + 0.209 ln(jetfuel) + 0.033 DIST + 0.018 ln(jetfuel)*DIST
# Road transport: ln(f/v)_{ckt} = -8.18 + 0.234 ln(w/v) + 0.862 ln(oil) + 1.21 ln(DIST) - 0.373ln(oil)*ln(DIST)
# Rail transport: ln(f/v)_{ckt} = -4.37 + 0.54 ln(w/v) + 0.079 ln(oil) - 0.90 ln(DIST) - 0.224 ln(oil)*ln(DIST)
# from: "The Wage Effects of Offshoring: Evidence from Danish Matched Worker-Firm Data" Hummels et al. 2013
DF2[,c('fv')] <- NA
DF2 <- within(DF2, fv[DF2$transp_mode == 'Water'] <- -4.16
+ 0.393*log(wv_ratio[DF2$transp_mode == 'Water'])
+ 0.351*log(oil[DF2$transp_mode == 'Water'])
+ 0.027*dist[DF2$transp_mode == 'Water']
+ 0.0063*log(oil[DF2$transp_mode == 'Water'])*(dist[DF2$transp_mode == 'Water']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Air'] <-  -3.80
+ 0.435*log(wv_ratio[DF2$transp_mode == 'Air'])
+ 0.209*log(jet[DF2$transp_mode == 'Air'])
+ 0.033*log(dist[DF2$transp_mode == 'Air'])
+ 0.018*log(jet[DF2$transp_mode == 'Air'])*(dist[DF2$transp_mode == 'Air']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Road'] <- -8.18
+ 0.234*log(wv_ratio[DF2$transp_mode == 'Road'])
+ 0.862*log(oil[DF2$transp_mode == 'Road'])
+ 1.21*log(dist[DF2$transp_mode == 'Road'])
- 0.373*log(oil[DF2$transp_mode == 'Road'])*log(dist[DF2$transp_mode == 'Road']))
DF2 <- within(DF2, fv[DF2$transp_mode == 'Rail'] <- -4.37
+ 0.54*log(wv_ratio[DF2$transp_mode == 'Rail'])
+ 0.079*log(oil[DF2$transp_mode == 'Rail'])
- 0.90*log(dist[DF2$transp_mode == 'Rail'])
- 0.224*log(oil[DF2$transp_mode == 'Rail'])*log(dist[DF2$transp_mode == 'Rail']))
# Getting Tau
DF2[,c('tau')] <- exp(DF2$fv)
#ad-valorem transport margins are CIF/FOB. If origin price is p and delivered price
#including margins is p*, and the advalorem margin is tau, then CIF/FOB gives you p*/p =
#(1+tau). CIF/FOB is in the column called "cif_fob"
# Getting the missing value, cif or fob. If the reported flow is cif, we need to find fob; if the reported flow is cif, we need
# to obtain cif. Therefore, if the reported values is cif, value2 will be the fob value and the other way around.
# cif = fob*(1+tau)
# fob = cif/(1+tau)
DF2[,c('value2')] <- NA
DF2 <- within(DF2, value2[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif']/(1 + tau[DF2$flow == 'cif']))
DF2 <- within(DF2, value2[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob']*(1 + tau[DF2$flow == 'fob']))
# We define a margin as cif - fob
DF2[,c('margin')] <- NA
DF2 <- within(DF2, margin[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'] - value2[DF2$flow == 'cif'])
DF2 <- within(DF2, margin[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'] - value_usd[DF2$flow == 'fob'])
# Getting the margin shares: margin/[cif_{total} - fob_{total}]
DF2[,c('cif_total')] <- 0
DF2[,c('fob_total')] <- 0
DF2 <- within(DF2, cif_total[DF2$flow == 'cif'] <- value_usd[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'cif'] <- value2[DF2$flow == 'cif'])
DF2 <- within(DF2, fob_total[DF2$flow == 'fob'] <- value_usd[DF2$flow == 'fob'])
DF2 <- within(DF2, cif_total[DF2$flow == 'fob'] <- value2[DF2$flow == 'fob'])
DF2[,c('odg')] <- paste(DF2$ISO_orig, DF2$ISO_dest, DF2$GTAP, sep="_")
DF2$cif_total <- ave(DF2$cif_total, DF2$odg, FUN=sum)
DF2$fob_total <- ave(DF2$fob_total, DF2$odg, FUN=sum)
DF2$share <- DF2$margin/(DF2$cif_total - DF2$fob_total)
#NANS <- DF2[is.na(DF2$share),]
DF2 <- DF2[!(is.na(DF2$share)),]
DF2 <- DF2[order(DF2$odg),]
row.names(DF2) <- NULL
# Column DF2$share is the one containing our margin shares. This one includes FOUR modes (air, sea, rail and road)
DT <- subset(DF2, select=c("ISO_orig", "ISO_dest", "value_usd", "weight_kg", "dist","GTAP", "transp_mode", "share", "wv_ratio"))
# We will repoort shares for air, water, and other (rail+ road)
DT[,c('mode')] <- NA
DT <- within(DT, mode[DT$transp_mode == "Rail"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Road"] <- 'Other')
DT <- within(DT, mode[DT$transp_mode == "Water"] <- 'Water')
DT <- within(DT, mode[DT$transp_mode == "Air"] <- 'Air')
DT$transp_mode <- NULL
DT[,c('odgm')] <- paste(DT$ISO_orig, DT$ISO_dest, DT$GTAP, DT$mode, sep="_")
DT$share <- ave(DT$share, DT$odgm, FUN=sum)
# "Other" is duplicated because it appears once where the mode was "rail" and once more where the mode was "road".
DT <- DT[!duplicated(DT$odgm), ]
DT <- DT[order(DT$odg),]
DT <- plyr::rename(DT, c("ISO_orig"="o", "ISO_dest"="d", "GTAP"="g", "mode"="m", "value_usd"="value_usd", "weight_kg"="weight_kg", "dist"="dist",
"share"="share", "wv_ratio"="wv_ratio", "odgm"="odgm" ))
write.csv(DT, "../DATA/output/OURSHARES.csv")
rm(list=ls())
workDir <- "C:/Users/jgnun/Documents/Purdue/Projects/Modal Shares/JGEA/R&R Sep 2018/Sumbited/codeForEstimation"
setwd(workDir)
########################################################################################
## Load R packages                                                                      ##
########################################################################################
library(xtable)
library(gdata)
library(dplyr)
library(reshape2)
library(tidyr)
library(ggplot2)
library(gridExtra)
library(bindrcpp)
############################
## GTAP bilateral margins and trade  (V 9.0)
vtwr.9.0 <- read.csv("../DATA/RawData/other/VTWR_9aY11.CSV")## sol( solution.dir = "../DATA/GTAP_Data/", solution.out = "VTWR_9aY11.HAR", csv.out = "VTWR_9aY11.CSV", header = "VTWR")
vxmd.9.0 <- read.csv("../DATA/RawData/other/VXMD_9aY11.CSV")
## Our new estimated shares:
nuvi <- read.csv("../DATA/output/OURSHARES.csv")
gtap.names <- read.csv("../DATA/RawData/other/gtapCodesComplete.csv")
bildist <- read.csv("../DATA/RawData/other/dist_cepii.csv") %>%
select( iso_o, iso_d, dist) %>%
dplyr::rename( o = iso_o, d = iso_d ) %>%
mutate( o = tolower(o),  d = tolower(d))
########################################################################################
## Clean and merge datasets                                                           ##
########################################################################################
vtwr <- rename.vars( data = vtwr.9.0,
from = c("MARG_COMM", "TRAD_COMM", "REG", "REG.1", "Value"),
to = c( "m", "g", "o", "d", "vtwr") )
vxmd <- rename.vars( data = vxmd.9.0,
from = c("TRAD_COMM", "REG", "REG.1", "Value"),
to = c("g", "o", "d", "vxmd") )
nuvi.mod <- nuvi %>%
mutate( o = tolower(o),
d = tolower(d),
m = recode( m , "Water" = "wtp", "Other" = "otp", "Air" =  "atp"),
value_MillUS = value_usd/1e6,
weight_ton = weight_kg/1e6) %>%
left_join( . , gtap.names[, c("g", "Code")],
by = ( "g") ) %>%
select( o, d, Code, m, value_MillUS, weight_ton, share, wv_ratio ) %>%
dplyr::rename( g = Code, wv_ratio.nuvi = wv_ratio, modal.share.nuvi = share )
## GTAP sectors to exclude: services, but also sugar cane and raw
## milk, where trade is marginal:
exclude <- c("c_b", "rmk", "ely","gdt","wtr","cns","trd","otp","wtp","atp","cmn","ofi","isr","obs","ros","osg","dwe")
## Weight-to-value ratios:
wv.ratios.by.g <- nuvi.mod %>%
group_by(o, d, g) %>%
summarise(wv.ratio = mean(wv_ratio.nuvi))
origins <- unique(vxmd$o) ## 141
destinations <- unique(vxmd$d) ## 141
sectors <- setdiff( unique(vxmd$g), exclude) ## 42 sectors
modes <- unique(vtwr$m) ## 3 modes
grid <- expand.grid( origins, destinations, sectors, modes) %>%
dplyr::rename( o = Var1, d = Var2, g = Var3, m = Var4 ) %>%
filter( o != d ) %>%
arrange(o , d, g, m)
## Add trade (the NAS here is because there is no trade: set them to
## zero):
dat <- left_join(grid, vxmd, by = c("o", "d", "g") ) %>%
left_join(. , vtwr, by = c("o", "d", "g", "m") )
dat[is.na(dat)] <- 0
## To calculate modal shares we need to add margins across modes,
## these are then merged back into dat, and the ratios (or modal
## shares) are taken:
transport.margins <- dat %>%
group_by(o, d, g)  %>%
summarise(transport.margin.gtap = sum(vtwr, na.rm = TRUE))
## Calculate modal shares:
dat <- left_join(dat, transport.margins, by = c("o", "d", "g") ) %>%
mutate(modal.share.gtap = vtwr/transport.margin.gtap)
