Numerical example for MLE for linear regression model I need to calculate the log-likelihood for a linear regression model in MATLAB. 
Although the theoretical result is well know and given in several sources, I want to find a numerical example so that I can check my code is correct.
Can anyone point me to one?
I realize that the parameters are the same as OLS (at least asymptotically) but its the actual log-likelihood I need.
 A: The generative model under OLS is that of: $y \sim N(X\beta, \sigma^2 I)$. Fitting the model is finding the $\beta$ and $\sigma$ that maximize your log-likelihood.
Assuming your general covariance is $K$ (here $K = \sigma^2 I$) the loglikehood $L$ is equal to :
$L = -\frac{N}{2}\log(2\pi) - \frac{1}{2} \log(|K|) - \frac{1}{2} (y-\hat{y})^T K^{-1}(y-\hat{y})$, where as I said above $N$ is your number of readings and now $\hat{y}$ are the model fitted values; $(y-\hat{y})$ are your model's residuals and $|K|$ denotes the determinant of the covariance matrix $K$.
Luckily for all of us the OLS log likelihood can also be expressed as :
$L = -\frac{N}{2}\log(2\pi) - N\log(\sigma) - \frac{1}{2\sigma^2} \sum (y-\hat{y})^2$. This and the above expression are equivalent. A bit of linear algebra can convince you for that.
OK, enough talk. Here is your numerical example in MATLAB:
clc; clear;                %just clear stuff

X = (1:.02:20)'; 
Y = cos(X);                %dependant variable
N = length(Y);             %number of readings
X_matrix = [ ones(N,1) X]; %make design matrix

[b] =  (X_matrix\ Y);      %solve the system
%[b2] = lscov(X_matrix,Y); %equivalent statement
Fitted = X_matrix * b;     %find fitted values
Residuals = Y - Fitted;    %find residuals
sigma = std(Residuals);    %find std.dev. of residuals
format long                %set it long for visual inspec. 
                           %Calculate L using simple expression
L_simple = -N*.5*log(2*pi) - N*log(sigma) - (1/(2*sigma.^2))*sum( Residuals.^2);
                           %set K covariance matrix
K_matrix =  eye(N) * sigma^2;
                           %Calculate L using generic expression
L_generic =-N*.5*log(2*pi) - sum(log(diag(chol(K_matrix)))) ...
           - .5*Residuals' / (K_matrix)* Residuals;
% ( sum(log(diag(chol(K_matrix)))) equals -.5*log(det(K)) 
% but it is far more stable.
L_generic 
% ans = -9.929305263221722e+02
L_generic - L_simple 
% ans =  3.410605131648481e-12 %cool they are practically the same

But hey, can you check this in R that we are quite sure it is works?
X = seq(1,20, by=.02)
Y = cos(X)
lm_test = lm(Y ~ X)
sigma = sd( residuals( lm_test))
logLik(lm_test)
#'log Lik.' -992.9303 (df=3)
#Check difference with MATLAB answer:
-9.929305263221722e+02 -  as.numeric(logLik(lm_test))
#[1] -0.0002630656 #very small difference, mostly due to numerics(*).
#Use the simple formula:
-N*.5*log(2*pi) - N*log(sigma) - (1/(2*sigma^2))*sum(  residuals( lm_test)^2)
#[1] -992.9305  #Cool it works as expected.

A standard reference free reference for all this would be the book Elements of Statistical Learning: Data Mining, Inference, and Prediction by Hastie, Tibshirani and Friedman. (Chapter 2.6)
(*) There are some small issues also about how you calculate the degrees of freedom for the variance; one uses (N-m) instead of just (N-1), m being the number of columns of your X matrix but you don't need to worry much about it at this point.
