# Calculate the coefficient of variation of production

# Load packages
library(rms)
library(ggplot2)
library(plyr)
library(reshape2)
library(xtable)

# Load World Production and Stocks variation data
Production <- read.csv("../Data/Data_Downloads/Production.csv",sep=";",header=T)
Production <- rename(Production,c("X"="Year"))
Production <- melt(Production,id.vars = "Year",variable.name = "Commodity")
Production$variable <- "Production"

StockVar <- read.csv("../Data/Data_Downloads/StockVar1000-1961-2011.csv",sep=";",header=T,
                     colClasses = c("integer",rep("numeric",10)))
StockVar <- melt(StockVar,id.vars = "Year",variable.name = "Commodity")
StockVar$variable <- "StockVar"

data <- rbind(Production,StockVar)
data <- dcast(data,Year+Commodity~variable)
# Calculate Consumption = Prod+StockVar
data$Consumption <- data$Production+data$StockVar
data <- melt(data,id.vars = c("Year","Commodity"),na.rm = TRUE)

# Copper consumption
CopperConsumption <- read.csv("../Data/Data_Downloads/CopperConsumption-1960-2013.csv",sep=";",header=T,
                     colClasses = c("integer","numeric"))
CopperConsumption$variable <- "Consumption"
CopperConsumption$Commodity <- "Copper_smelter"
CopperConsumption <- rename(CopperConsumption,c("Copper"="value"))

data <- rbind(data,CopperConsumption)
levels(data$Commodity)[4] <- "Copper"
levels(data$Commodity)[8] <- "Palm oil"
levels(data$Commodity)[12] <- "Tin"
data <- subset(data,Year>=1961 & Year<=2011)
data <- droplevels(data)

# Choose knots
knots3 <- c(1966,1986,2006)
knots4 <- c(1964,1979,1994,2008)
knots5 <- c(1964,1975,1986,1998,2008)

# Plot all quantities
AllQuantities <- ggplot(subset(data,variable!="StockVar"),
                        aes(x = Year, y = value/1E6, colour = variable)) +
                        geom_line(size = 1) +
                        facet_wrap(~Commodity,scales = "free_y",ncol=3) +
                        theme_bw() +
                        ylab("Quantities (million ton)") +
                        guides(colour = guide_legend(title = NULL), linetype = guide_legend(title = NULL)) +
                        theme(legend.position = "bottom", legend.key = element_blank())

ggsave(filename = "../Figures/Times/AllQuantities.pdf",
       plot = AllQuantities,
       units = "in",
       width = 12,
       height = 15,
       pointsize = 9,
       family = "Times")

# Detrend the quantities using the rcs and calculate CV
CalculateCV <- function(subdata) {
    subdata <- subset(subdata,variable=="Production")
    subdata.rcs <- lm(log(value)~ns(Year,knots=knots5[2:4],Boundary.knots = knots5[c(1,5)]),
                      subdata)
    subdata <- merge(subdata,residuals(subdata.rcs),by = 0) # Merge by row names
    ddply(subdata,.(variable),summarize, CV = sd(y))
}
Stat <- dcast(ddply(data,.(Commodity),CalculateCV),Commodity~variable)

# Plot the trend
logdata <- transform(subset(data,variable!="StockVar"),
                     value = log(value))
logdata <- ddply(logdata,.(Commodity),mutate,
                 trend = fitted(lm(value~ns(Year,knots=knots5[2:4],Boundary.knots = knots5[c(1,5)]))))
ggplot(subset(logdata,variable!="StockVar"),
       aes(x = Year, y = value, colour = variable)) +
       geom_line(size = 1) +
       geom_line(aes(x=Year,y=trend),size = 0.5,colour = "black")+
       facet_wrap(~Commodity,scales = "free_y",ncol=3) +
        theme_bw() +
        ylab("Quantities (log of million ton)") +
        guides(colour = guide_legend(title = NULL), linetype = guide_legend(title = NULL)) +
        theme(legend.position = "bottom", legend.key = element_blank())

# Save the CV
write.table(Stat, file = "../Results/CVQuantities.csv", quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = TRUE, qmethod = c("escape", "double"),
            fileEncoding = "")

# Export to latex
StatTex <- Stat
StatTex$Production <- formatC(StatTex$Production*100,digits = 2,flag = "#",format = "f")
StatTex$Commodity <- as.character(StatTex$Commodity)
StatTex$Consumption <- NULL
(StatTex <- cbind(StatTex[1:7,],rbind(StatTex[8:13,],c("","",""))))

print(xtable(StatTex),
      file = "../Article/Tables/SDProdConso.tex",
      only.contents = TRUE,
      include.rownames = FALSE,
      include.colnames = FALSE,
      hline.after = NULL)
