car = read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/h14/data/car.txt", sep=",", header=T)
head(car)

##############################################################################
####Exercise 3.2c) from the book
##############################################################################

#We study the dataset car
claimSizeAll = as.numeric(car$claimcst0)                 # A numeric vector, between 0 and infinity
claimSize    = claimSizeAll[(claimSizeAll > 0)]/1000     # Exclude(!) all data where claim_i=0, scale by 1/1000 (as explained in the book on page 15 in the book)

#let claimSize be gamma(nu, mu) as equation 2.6 in the Jong and Heller
#find the MLE functions:
MLE = function(params){
            nu=params[1]
            mu=params[2]
            n=length(claimSize)
            LogLik=(nu-1)*sum(log(claimSize)) - n*log(gamma(nu)) + n*nu*log(nu/mu) - nu/mu*sum(claimSize)
            return(-LogLik)
}

init = c(2,2)
est = optim(par=init, fn=MLE)$par      #\hat{nu}, \hat{mu}
nu = est[1]
mu = est[2]
#optim is a minimiser: minimize -LogLik i.e. maximizes LogLik.


#re-parameterise to alpha=shape, beta=rate  as explained in the book on page 28, because this is the parameterisation R uses.
alpha_est =  nu
beta_est  =  nu/mu

#PLOT
par(mfrow=c(1,2))
hist(claimSize, xlim=c(0,15), 100)        #For clarity the horizontal axis is truncated at 15.
x = seq(from=1, to=15, by=0.1)
plot(x, dgamma(x, shape=alpha_est, rate = beta_est), ylab="P(claimSize)", col="red", type="l")




##############################################################################
####Exercise 3.3d) from the book    #let claimSize be IG(mu, sigma^2)
##############################################################################

#find MLE:
MLE = function(params){
            mu=params[1]
            sigma2=params[2]
            LogLik=sum(-log(sqrt(2*pi*claimSize^3*sigma2)) - 1/(2*claimSize)*((claimSize-mu)/(mu*sigma2))^2)
        return(-LogLik)
}

init = c(2,2)
est = optim(init, MLE)$par      #\hat{nu}, \hat{mu}
mu_est = est[1]
sigma2_est = est[2]

y = seq(from=1, to=15, by=0.5)
fy = 1/sqrt(2*pi*y^3*sigma2_est)*exp(-1/(2*y)*((y-mu_est)/(mu_est*sigma2_est))^2)        #calculate the IG pdf.
lines(y, fy, ylab="P(claimSize)", col="blue", type="l")

legend("topleft",legend=c("gamma","IG"),col=c("red","blue"),lty=1)