# Plot the price function for the model estimated without trend for
# wheat along with the observed prices and the ergodic distribution of
# availability

source("functions.R")
library(R.matlab)

Wheat <- readMat("../Results/Wheat.mat")

PlotWheatPrice <- function() {
    par(fig = c(0,1,0,0.77),
        new=TRUE,
        las = 1,
        mar = c(4.1,4.1,0,0.5))
    plot(Wheat$s,
         Wheat$x[,2],
         type = 'l',
         xlab = 'Availability',
         ylab = 'Price',
         yaxs = "i",
         axes = FALSE,
         xlim = c(-0.5,max(Wheat$Aobs)),
         ylim = c(0,5))
    points(Wheat$Aobs,
           Wheat$Pobs,
           pch = '+',
           col = "blue")
    abline(h = mean(Wheat$Pobs),
           lty = 2)
    abline(h = Wheat$a,
           lty = 3)
    abline(h = Wheat$pstar,
           lty = 4)
    text(x = 12, y = mean(Wheat$Pobs), labels = "Sample average",adj = c(1,-0.3))
    text(x = 12, y = Wheat$a, labels = "Ergodic mean",adj = c(1,-0.3))
    text(x = 12, y = Wheat$pstar, labels = expression(italic(P^"*")),adj = c(1,-0.3))
    axis(side = 1, at=seq(-2,16,by=2))
    axis(side = 2, at=seq(-1,6,by=1))

    par(fig = c(0,1,0.77,1),
        new=TRUE,
        mar = c(0,4.1,0.2,0.5))
    hist(Wheat$Asim,
         125,
         freq=FALSE,
         xlim = c(-0.5,max(Wheat$Aobs)),
         main = "",
         axes = FALSE,
         ylab = "",
         xlab = "",
         border = "#EEE8CD")
    lines(density(Wheat$Asim))
}
PlotWheatPrice()

FiguresCreation(file = c("../Figures/","WheatPrice"),
                quote(PlotWheatPrice),
                width=8,height=6,pointsize=11)
