I am trying to implement a multivariate linear regression model by ML estimation, however I ran into some discrepancies.

An assumption of the model is that the residual is a standard normal variable.

I followed this(page 6) approach, but the procedure should be clear anyway:

$$ \text{pdf}_X (x,\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$

$$ \text{pdf}_X (x_1,...,x_n,\mu,\sigma) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x_i-\mu)^2}{2\sigma^2}} = (\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}} $$

Taking the log gives us the log likelihood:

$$ log((\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}}) = -\frac{n}{2}log(2\pi) -\frac{n}{2}log(\sigma^2) - \frac{\sum_{i=1}^n(y - X\beta)^2}{2\sigma^2} $$

Translating this into R Code:

set.seed(1)

# generate random data 
X <- cbind(rep(1, 100/4), matrix(rnorm(100), ncol=4, byrow=T))
#true beta 
beta <- c(25, 1,3,5,7)
# linear model
y <- X%*%beta + rnorm(25)

llike <- function(par, y, X){
  m <- nrow(X)
  n <- ncol(X)
  beta <- par[1:n]
  sigma_sq <- par[1+n]
  e <- y - X%*%beta
  loglike <- -(m/2)*log(2*pi) - (m/2)*log(sigma_sq) - ((t(e)%*%e) / (2*sigma_sq))
  return(-loglike)
}

res <- optim(par=c(rep(1, 6)), llike, method="BFGS", hessian=T, y=y, X=X)

res$par
24.189 1.209 3.729 5.489 7.102 8788.176

The last term is the variance it seems excessive, but keep in mind that the "rnorm" function generates normal distributed data, with mean = 0 and standard deviation = 1. Furthermore we have 25 of those so the estimated standard deviation would actually be:

 sqrt(res$par[6])/25
 3.749

Now when I calculate the logliklihood:

-llike(c(res$par[1:5], sqrt(res$par[6])/25), y, X)
-45.56091

compared to:

logLik(lm(y ~ X-1))
-42.238

So there are two issues that I have:

1. The estimates between the optim and lm are different
2. the loglikelihood is differes.

Am I missing something with the implementation? I really need to know if there is a mistake, since this model is the bases for Lasso and Ridge Regression which are supposed to follow.

I am trying to implement a multivariate linear regression model by ML estimation, however I ran into some discrepancies.

An assumption of the model is that the residual is a standard normal variable.

I followed this(page 6) approach, but the procedure should be clear anyway:

$$ \text{pdf}_X (x,\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$

$$ \text{pdf}_X (x_1,...,x_n,\mu,\sigma) = \prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x_i-\mu)^2}{2\sigma^2}} = (\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}} $$

Taking the log gives us the log likelihood:

$$ log((\frac{1}{\sqrt{2\pi\sigma^2}})^ne^{-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}}) = -\frac{n}{2}log(2\pi) -\frac{n}{2}log(\sigma^2) - \frac{\sum_{i=1}^n(y - X\beta)^2}{2\sigma^2} $$

Translating this into R Code:

set.seed(1)

# generate random data 
X <- cbind(rep(1, 100/4), matrix(rnorm(100), ncol=4, byrow=T))
#true beta 
beta <- c(25, 1,3,5,7)
# linear model
y <- X%*%beta + rnorm(25)

llike <- function(par, y, X){
  m <- nrow(X)
  n <- ncol(X)
  beta <- par[1:n]
  sigma_sq <- par[1+n]
  e <- y - X%*%beta
  loglike <- -(m/2)*log(2*pi) - (m/2)*log(sigma_sq) - ((t(e)%*%e) / (2*sigma_sq))
  return(-loglike)
}

res <- optim(par=c(rep(1, 6)), llike, method="BFGS", hessian=T, y=y, X=X)

res$par
24.189 1.209 3.729 5.489 7.102 8788.176

The last term is the variance it seems excessive, but keep in mind that the "rnorm" function generates normal distributed data, with mean = 0 and standard deviation = 1. Furthermore we have 25 of those so the estimated standard deviation would actually be:

 sqrt(res$par[6])/25
 3.749

Now when I calculate the logliklihood:

-llike(c(res$par[1:5], sqrt(res$par[6])/25), y, X)
-45.56091

compared to:

logLik(lm(y ~ X-1))
-42.238

So there are two issues that I have:

1. The estimates between the optim and lm are different
2. the loglikelihood is differes.

Am I missing something with the implementation? I really need to know if there is a mistake, since this model is the bases for Lasso and Ridge Regression which are supposed to follow.

EDIT

Thanks the suggestions provided by Tim and the solution from here it tunes out there were two mistakes:

1. in the llike function

   sigma_sq <- par[1+n] #replace by sigma_sq <- par[1+n]^2

2. the optimization algorithm:

Change from "BFGS" to "L-BFGS-B" and obtain the desired results.