### Using maxLik package in R
### By: Mauricio Sarrias


## ----loadmaxLik-----------------------------------------------------------------------------------------------------------
library("maxLik")


## ----llog1----------------------------------------------------------------------------------------------------------------
linear.ml <- function(param, y, X){
  n <- nrow(X)
  k <- ncol(X)
  beta <- param[1:k] 
  sig.sq  <- param[k + 1]
  resi  <- y - tcrossprod(X, t(beta))
  LL    <- -0.5 * n * log(2 * pi * sig.sq) - 0.5 * crossprod(resi) / sig.sq 
  return(LL)
}


## ----art_data-------------------------------------------------------------------------------------------------------------
# Create artifitial database
set.seed(123)
N  <- 5000     # number of individuals
x1 <- runif(N)
x2 <- rnorm(N)
y <- 1 + 1 * x1 + 1 * x2 + rnorm(N, mean = 0, sd = 2) # e~N(0, 2^2)
X <- cbind(1, x1, x2) # put data in matrix form


## ----initial-values-------------------------------------------------------------------------------------------------------
start <- c(0, 0, 0, 1)
names(start) <- c("constant", "x1", "x2", "sigma2")


## ----out1-----------------------------------------------------------------------------------------------------------------
# Estimate the model by BFGS
out1 <- maxLik(logLik = linear.ml, start = start, 
               y = y, X = X, method = 'bfgs')


## ----out1sum--------------------------------------------------------------------------------------------------------------
summary(out1)


## ----coef-----------------------------------------------------------------------------------------------------------------
coef(out1)
stdEr(out1)


## ----avcov----------------------------------------------------------------------------------------------------------------
avov <- -solve(out1$hessian) #extract the Hessian at theta hat
se   <- sqrt(diag(avov))
se
all.equal(se, stdEr(out1))


## ----mlgr-----------------------------------------------------------------------------------------------------------------
linear.mlb <- function(param, y, X){
  k <- ncol(X)
  beta <- param[1:k] 
  sig.sq  <- param[k + 1]
  resi  <- y - tcrossprod(X, t(beta))
  # Log-likelihood function for each individual
  Li    <-   -0.5 * log(2 * pi * sig.sq) - 0.5 * resi ^ 2 / sig.sq 
  
  #Score function for each individual
  grad.beta <- (1 / sig.sq) * X * drop(resi)
  grad.sig  <- -1/(2 * sig.sq) + resi ^ 2 / (2 * sig.sq ^ 2)
  grad <- cbind(grad.beta, grad.sig)
  attr(Li, 'gradient') <- grad
  Li
}


## ----sout2----------------------------------------------------------------------------------------------------------------
out2 <- maxLik(logLik = linear.mlb, start = start, 
               y = y, X = X, method = 'bhhh')
summary(out2)


## ----check----------------------------------------------------------------------------------------------------------------
all.equal(coef(out1), coef(out2))
all.equal(stdEr(out1), stdEr(out2))


## ----sebhhh---------------------------------------------------------------------------------------------------------------
G <- out2$gradientObs
GG <- crossprod(G)
se_bhhh <- sqrt(diag(solve(GG)))
se_bhhh
all.equal(se_bhhh, stdEr(out2))


## ----code3, size = 'scriptsize'-------------------------------------------------------------------------------------------
linear.mlc <- function(param, y, X){
  k <- ncol(X)
  beta <- param[1:k] 
  sig.sq  <- param[k + 1]
  if (sig.sq <= 0) return(NA) # This signals to the optimazer
  # that the attempted parameter value was out of range, and forces
  # it to find a new one closer to the previous values
  resi  <- y - tcrossprod(X, t(beta))
  Li    <-   -0.5 * log(2 * pi * sig.sq) - 0.5 * resi ^ 2 / sig.sq 

  #Gradient
  grad.beta <- (1 / sig.sq) * X * drop(resi)
  grad.sig  <- -1/(2 * sig.sq) + resi ^ 2 / (2 * sig.sq ^ 2)
  grad <- cbind(grad.beta, grad.sig)
  attr(Li, 'gradient') <- grad

  
  #Hessian
  H <- matrix(NA, nrow = k + 1, ncol = k + 1)
  H[1:k, 1:k] <- -(1 / (sig.sq)) * crossprod(X)
  H[k + 1, k + 1] <- 1 / (2 * sig.sq ^ 2) - (1 / (sig.sq ^ 3)) * crossprod(resi)
  H[1:k, k + 1] <- -(1 / (sig.sq ^ 2))  * crossprod(X, resi)
  H[k + 1, 1:k] <- -(1 / (sig.sq ^ 2))  * crossprod(resi, X)
  attr(Li, 'hessian') <- H
  Li
}



## ----nr-------------------------------------------------------------------------------------------------------------------
out3 <- maxLik(logLik = linear.mlc, start = start, 
               y = y, X = X, method = 'nr')
summary(out3)


## ----ols, size = 'scriptsize'---------------------------------------------------------------------------------------------
colnames(X) <- c("constant", "x1", "x2")
ols <- lm(y ~ X - 1)
summary(ols)


## -------------------------------------------------------------------------------------------------------------------------
delta <- seq(-1, 4, .01)
grad.values <- as.numeric(lapply(delta, function(x) {
                            colSums(attr(linear.mlb(param <- c(1, 1, x, 4), 
                                                    y = y, 
                                                    X = X), 
                                         'gradient')
                                    )[[3]]     
                                      }
                                )
                          )

logl.values <- as.numeric(lapply(delta, function(x) {
   sum(linear.mlb(param <- c(1, 1, x, 4),y = y,X = X))     
}))


## ----plotLL1, echo = FALSE, out.width='3in'-------------------------------------------------------------------------------
#par(mar = c(5, 4, 4, 5) + .1)
plot(delta, logl.values, 
     type = 'l', 
     lty  = 'solid',
     lwd  = 3, 
     col  = 'red', 
     xlab = expression(beta[2]), 
     ylab = 'Log Likelihood')
par(new = TRUE)
plot(delta, grad.values, 
     type  = 'l',
     lty   = 'dashed',
     lwd   = 3, 
     col   = 'blue', 
     xlab  = "", 
     ylab  = "",
     xaxt  = "n",
     yaxt = "n")
axis(4)
mtext("Gradient", side = 4, line = 3)
legend("bottomleft",
       col = c("red","blue"),
       lty = c("solid", "dashed"),
       legend = c("Log-Likelihood","Gradient"),
       cex = 1.2, 
       #y.intersp = 0.2, 
       bty = "o",
       lwd = 2)


## ----plotLL2, echo = FALSE, out.width='3in'-------------------------------------------------------------------------------
b_1 <- seq(-2, 4, .1)
b_2 <- seq(-2, 4, .1)

#grid of values
g <-   expand.grid(beta1 = b_1, beta = b_2)
g$z <- apply(g, 1, function(x) sum(linear.mlb(param <- c(1, x[1], x[2], 4),y = y,X = X)) )

# reorganize z as a matrix same shape as grid
zmat <- matrix(g$z, nrow = length(b_1))
contour(b_1, b_2, zmat,
        xlab = expression(beta[1]),
        ylab = expression(beta[2]))
points(1, 1, pch = 16, col = 2, cex = 1.2)
abline(h = 1)
abline(v = 1)


## ----plotLL3, echo = FALSE, out.width='3in'-------------------------------------------------------------------------------
persp(b_1, b_2, zmat, 
      xlab  = expression(beta[1]), 
      ylab  = expression(beta[2]), 
      zlab  = "log-likelihood", 
      theta = 30, 
      phi   = 30)