sandbox/analysis_prod_trade_stocks.R

In SWS-Methodology/faoswsStock: Package to perform the imputation of stock changes for the FAO food balance sheets

# This code analysis for production, trade and stocks by cereals, pulses and sugar.

## load the library
library(faosws)
library(data.table)
library(ggplot2)
library(faoswsUtil)

## set up for the test environment and parameters
# initialYear = 1961 # Change documentation if this changes!

R_SWS_SHARE_PATH <- Sys.getenv("R_SWS_SHARE_PATH")

# This return FALSE if on the Statistical Working System
if(CheckDebug()){
  
  message("Not on server, so setting up environment...")
  
  library(faoswsModules)
  SETTINGS <- ReadSettings("modules/impute_stocks/sws.yml")
  
  # If you're not on the system, your settings will overwrite any others
  R_SWS_SHARE_PATH <- SETTINGS[["share"]]
  
  # Define where your certificates are stored
  SetClientFiles(SETTINGS[["certdir"]])
  
  # Get session information from SWS. Token must be obtained from web interface
  GetTestEnvironment(baseUrl = SETTINGS[["server"]],
                     token = SETTINGS[["token"]])
  
}

elements = c("51", "61", "71", "91")

# 51 production (Mt)
# 61 imports (Mt)
# 71 Variation Intial Exstenc (Mt) 
# 91 exports (Mt)  
# 101 feed (Mt)
# 131 processed (Mt)
# 141 food (Mt)

# Data
key = DatasetKey(domain = "faostat_one", dataset = "FS1_SUA_UPD", dimensions = list(
  Dimension(name = "geographicAreaFS", keys = GetCodeList("faostat_one", "FS1_SUA_UPD", "geographicAreaFS")[, code]),
  Dimension(name = "measuredElementFS", keys = elements),
  Dimension(name = "timePointYears", keys = as.character(2005:2012)),
  Dimension(name = "measuredItemFS", keys = GetCodeList("faostat_one", "FS1_SUA_UPD", "measuredItemFS")[, code]))
)
data = GetData(key)
data[, geographicAreaM49 := fs2m49(as.character(geographicAreaFS))]
data[, measuredItemCPC := fcl2cpc(as.character(formatC(as.numeric(measuredItemFS), width = 4,
                                                       flag = "0")))]
# data = data[!is.na(geographicAreaM49) & 
#               !(flagFaostat %in% c("M", "T"))]

data[, c("geographicAreaFS", "measuredItemFS", "flagFaostat") := NULL]
setcolorder(data, c("geographicAreaM49", "measuredItemCPC", "measuredElementFS", 
                    "timePointYears", "Value"))

# data[grepl("Hen", measuredItemCPC_description)]

data = data[, list(Value = sum(Value)), 
            by = list(geographicAreaM49, measuredItemCPC, 
                      measuredElementFS, timePointYears)]

item <- ReadDatatable("fbs_stocks_comm_codes")
item <- item[fbs_key == "measuredItemCPC", fbs_code]

# Select just cereals, pulses and refined sugar

data <- data[measuredItemCPC %in% item]
data[substr(measuredItemCPC, 1, 3) == "011", group := "cereals"]
data[substr(measuredItemCPC, 1, 3) == "017", group := "pulses"]
data[substr(measuredItemCPC, 1, 3) == "235", group := "refined_sugar"]
data

## More data manipulation 
data <- dcast.data.table(data, geographicAreaM49 + measuredItemCPC + timePointYears + group ~ measuredElementFS,
                         value.var = "Value")

setnames(data, old=c("51", "61", "71", "91"),
         new=c("production", "imports", "variationIntialExstenc", "exports"))

## Checking

for(cname in c("production", "imports", "variationIntialExstenc", "exports")){
  data[is.na(get(cname)), c(cname) := 0]
}

data[, netTrade := imports - exports]
data

## Country groups
dlpath <- file.path("C:", "Users", "caetano", "Documents", "Github", "faoswsFeed",
                    "data-raw", "IR_factor", "class.csv")

countryGroup <- fread(dlpath)
countryGroup[, .N, GroupName]

countryGroup[grepl("income", GroupName), .N, GroupName]

countryIncomeGroup <- countryGroup[GroupCode %in% c("HIC", "LIC", "UMC", "LMC"), ]
# countryRegionGroup <- countryGroup[GroupCode %in% c("NAC", "LCN", "AFR", "EAS", "EUU"),]

# library(countrycode)
# ?countrycode

countryIncomeGroup[, geographicAreaM49 := as.character(countrycode(CountryCode, "wb", "iso3n"))]

# Sudan has the wrong name (it should be former Sudan)
countryIncomeGroup[geographicAreaM49 == "736", CountryName := "Sudan (former)"]
# China should be 1248
countryIncomeGroup[geographicAreaM49 == "156", geographicAreaM49 := "1248"]
#Exclude Channel Islands and Kosovo (not separately recognised by the UN)
countryIncomeGroup <- countryIncomeGroup[!is.na(geographicAreaM49)]
countryIncomeGroup

# Merge

data <- merge(data, countryIncomeGroup[, c("geographicAreaM49", "GroupName"), with = F], 
      by = "geographicAreaM49", all.x=T)

data <- data[!is.na(geographicAreaM49)]
data <- data[!is.na(GroupName)]

tabGroup <- data[, list(production = sum(production, na.rm = T),
            imports = sum(imports, na.rm = T),
            exports = sum(exports, na.rm = T),
            deltaStocks = sum(variationIntialExstenc, na.rm = T)),
     by=list(GroupName, group, timePointYears)]

tabGroup <- tabGroup[, trade := imports - exports]

tabGroup[timePointYears == 2010 & group == "refined_sugar"]
getwd()
# saving data for the r markdown file
# write.csv(tabGroup, "sandbox/tabGroup.csv", row.names = F)



tabGroup$GroupName <- factor(tabGroup$GroupName, 
                        levels = c("High income", "Upper middle income", 
                                   "Lower middle income", "Low income"))

ggplot(data=tabGroup,
       aes(x=as.numeric(timePointYears), y=production, group=GroupName, col = GroupName)) +
  geom_line(aes(), stat = "identity", position=position_dodge(), size=1.5) +
  facet_wrap(~ group
             , scales = "free"
  ) +
  theme(legend.title=element_blank(), strip.text = element_text(size=7)) +
  ylab('Year') + xlab('Production') +
  scale_x_continuous(lim=c(2005, 2012), breaks=seq(2005, 2012, 1)) +
  theme_fao(45)


ggplot(data=tabGroup,
       aes(x=as.numeric(timePointYears), y=trade, group=GroupName, col = GroupName)) +
  geom_line(aes(), stat = "identity", position=position_dodge(), size=1.5) +
  facet_wrap(~ group
             , scales = "free"
  ) +
  theme(legend.title=element_blank(), strip.text = element_text(size=7)) +
  ylab('Year') + xlab('Imports - Exports') +
  scale_x_continuous(lim=c(2005, 2012), breaks=seq(2005, 2012, 1)) +
  geom_hline(yintercept=0) +
  theme_fao(45)