# Plot all prices with their respective trends
source("functions.R")
library(R.matlab)
library(reshape2)
library(ggplot2)

# Data
## Import data from Matlab
Without <- readMat("../Results/WithoutPrices.mat")$Prices
Linear <- readMat("../Results/LinearPrices.mat")$Prices
Quadratic <- readMat("../Results/QuadraticPrices.mat")$Prices
Cubic <- readMat("../Results/CubicPrices.mat")$Prices
RCS3 <- readMat("../Results/RCS3Prices.mat")$Prices
RCS4 <- readMat("../Results/RCS4Prices.mat")$Prices

Year <- 1900:2011

i <- 8

PlotWithoutTrendPrices <- function() {
  par(mar = c(4.1,4.1,0.5,0.5),
      las = 1)
  ii <- seq(from = 1, by = 3, length.out = 13)[i]
  Prices <- Without[[ii+1]]
  matplot(Year,
          Prices,
          ylab = "Real prices",
          type = "l",
          xaxs = "i",
          lwd = 1.5,
          lty = 1,
          col = "black")
  abline(h=Without[[ii+2]],col = "lightblue",lwd = 1, lty = 2)
}
PlotWithoutTrendPrices()

PlotTrendPrices <- function(data = Linear,LOG = TRUE) {
  par(mar = c(4.1,4.1,0.5,0.5),
      las = 1)
  ii <- seq(from = 1, by = 5, length.out = 13)[i]
  Prices <- cbind(data[[ii+1]],data[[ii+2]],data[[ii+3]])
  if(LOG) {
      ylab <- "Prices (log)"
      y <- log(Prices)
  }
  else {
      ylab <- "Prices"
      y <- Prices
  }
  matplot(Year,
          y,
          type = "l",
          xaxs = "i",
          ylab = ylab,
          lwd = c(1.5,1.5,1),
          lty = c(1,1,3),
          col = c("black","blue","grey15"))
  if(LOG) abline(h=log(data[[ii+4]]),col = "lightblue",lwd = 1, lty = 2)
  else    abline(h=data[[ii+4]],col = "lightblue",lwd = 1, lty = 2)
  legend(x = "topright",
         bty = "n",
         lty = c(1,1,3,2),
         col = c("black","blue","grey15","lightblue"),
         legend = c("Real price","Detrended price","Trend",
                    expression(paste("Cutoff price ",group('(',italic(p^"*"),')')))))
}
PlotTrendPrices()

i <- 1
PlotTrendPrices(data = RCS3,LOG = FALSE)

PlotAllTrendRCS <- function() {
  par(mar = c(4.1,4.1,0.5,0.5),
      las = 1)
  ii <- seq(from = 1, by = 5, length.out = 13)[i]
  meanPobs <- mean(Linear[[ii+1]])
  Pobs <- Linear[[ii+1]]
  ToMinimize <- function(coef,Trend) sum((Pobs-coef*Trend)^2)
  coeflin <- optim(meanPobs/mean(Linear[[ii+3]]),function(coef) ToMinimize(coef,Linear[[ii+3]]))$par
  coefRCS3 <- optim(meanPobs/mean(RCS3[[ii+3]]),function(coef) ToMinimize(coef,RCS3[[ii+3]]))$par
  coefRCS4 <- optim(meanPobs/mean(RCS4[[ii+3]]),function(coef) ToMinimize(coef,RCS4[[ii+3]]))$par
  Prices <- cbind(Linear[[ii+1]],
                  Linear[[ii+3]]*coeflin,
                  RCS3[[ii+3]]*coefRCS3,
                  RCS4[[ii+3]]*coefRCS4)
  matplot(Year,
          Prices,
          type = "l",
          xaxs = "i",
          ylab = "Prices",
          lwd = 1.5,
          lty = 1:4,
          col = 1:4)
  legend(x = "topright",
         bty = "n",
         lty = 1:4,
         col = 1:4,
         lwd = 1.5,
         legend = c("Real price","Linear trend","3-knot spline trend","4-knot spline trend"))
}
PlotAllTrendRCS()

# Gather all prices in one data.frame to plot them with ggplot
for(i in 1:13) {
  ii <- seq(from = 1, by = 5, length.out = 13)[i]
  meanPobs <- mean(Linear[[ii+1]])
  Pobs <- Linear[[ii+1]]
  ToMinimize <- function(coef,Trend) sum((Pobs-coef*Trend)^2)
  coeflin <- optim(meanPobs/mean(Linear[[ii+3]]),function(coef) ToMinimize(coef,Linear[[ii+3]]))$par
  coefRCS3 <- optim(meanPobs/mean(RCS3[[ii+3]]),function(coef) ToMinimize(coef,RCS3[[ii+3]]))$par
  coefRCS4 <- optim(meanPobs/mean(RCS4[[ii+3]]),function(coef) ToMinimize(coef,RCS4[[ii+3]]))$par
  Prices <- cbind(Linear[[ii+1]],
                  Linear[[ii+3]]*coeflin,
                  RCS3[[ii+3]]*coefRCS3,
                  RCS4[[ii+3]]*coefRCS4)
  colnames(Prices) <- c("Without","Linear","RCS3","RCS4")
  Prices <- as.data.frame(Prices)
  Prices$Year <- Year
  Prices$Commodity <- Without[[(i-1)*3+1]][1][1]
  if(i==1) AllPrices <- Prices
  else AllPrices <- rbind(AllPrices,Prices)
}
AllPrices$Commodity <- as.factor(AllPrices$Commodity)
AllPrices <- melt(AllPrices,measure.vars = c("Without","Linear","RCS3","RCS4"),
                  variable.name = "Trend")

levels(AllPrices$Commodity)[8] <- "Palm oil"
levels(AllPrices$Trend) <- c("Real price","Linear trend","3-knot spline trend","4-knot spline trend")

PlotAllPrices <- function() {
    ggplot(AllPrices,aes(x = Year, y = value,colour = Trend)) +
        geom_line(size = 1,aes(linetype=Trend)) +
        facet_wrap(~Commodity,ncol=3,scales = "free_y") +
        theme_bw() +
        ylab("Prices") +
        scale_colour_manual(values = c("black","red","blue", "darkgreen")) +
        guides(colour = guide_legend(title = NULL), linetype = guide_legend(title = NULL)) +
        theme(legend.position = "bottom", legend.key = element_blank(), text = element_text(size = 20))
}
PlotAllPrices()

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

PlotAllTrendPoly <- function() {
  par(mar = c(4.1,4.1,0.5,0.5),
      las = 1)
  ii <- seq(from = 1, by = 5, length.out = 13)[i]
  meanPobs <- mean(Linear[[ii+1]])
  Pobs <- Linear[[ii+1]]
  ToMinimize <- function(coef,Trend) sum((Pobs-coef*Trend)^2)
  coeflin <- optim(meanPobs/mean(Linear[[ii+3]]),function(coef) ToMinimize(coef,Linear[[ii+3]]))$par
  coefQuadratic <- optim(meanPobs/mean(Quadratic[[ii+3]]),function(coef) ToMinimize(coef,Quadratic[[ii+3]]))$par
  coefCubic <- optim(meanPobs/mean(Cubic[[ii+3]]),function(coef) ToMinimize(coef,Cubic[[ii+3]]))$par
  Prices <- cbind(Linear[[ii+1]],
                  Linear[[ii+3]]*coeflin,
                  Quadratic[[ii+3]]*coefQuadratic,
                  Cubic[[ii+3]]*coefCubic)

  matplot(Year,
          Prices,
          type = "l",
          xaxs = "i",
          ylab = "Prices",
          lwd = 1.5,
          lty = 1:4,
          col = 1:4)
  legend(x = "topright",
         bty = "n",
         lty = 1:4,
         col = 1:4,
         lwd = 1.5,
         legend = c("Real price","Linear trend","Quadratic trend","Cubic trend"))
}
PlotAllTrendPoly()


for(i in 1:13) {
    FiguresCreation(file = c("../Figures/",
                    paste(c("Without/",Without[[(i-1)*3+1]][1][1],"Price"),collapse="")),
                    quote(PlotWithoutTrendPrices),
                    width = 4,
                    height = 3,
                    pointsize = 9)
    for(trend in c("Linear","Quadratic","Cubic","RCS3","RCS4")) {
        plotfun <- function() PlotTrendPrices(data = get(trend),LOG = FALSE)
        FiguresCreation(file = c("../Figures/",
                        paste(c(trend,"/",get(trend)[[(i-1)*5+1]][1][1],"Price"),collapse="")),
                        quote(plotfun),
                        width = 4,
                        height = 3,
                        pointsize = 9)
  }
    FiguresCreation(file = c("../Figures/",
                    paste(c("AllTrendRCS/",Without[[(i-1)*3+1]][1][1],"Price"),collapse="")),
                    quote(PlotAllTrendRCS),
                    width = 4,
                    height = 3,
                    pointsize = 9)
    FiguresCreation(file = c("../Figures/",
                    paste(c("AllTrendPoly/",Without[[(i-1)*3+1]][1][1],"Price"),collapse="")),
                    quote(PlotAllTrendPoly),
                    width = 4,
                    height = 3,
                    pointsize = 9)
}

