How to calculate the Lagrange multiplier test statistic?

#1
I am trying to do this on Matlab.

n = 200;

%The model has 4 variables X1, X2, X3, X4

X1=rand(n,1); X2=rand(n,1); X3=rand(n,1); X4=rand(n,1); X=[X1,X2,X3,X4];

beta1=-0.6; beta2=-0.5; beta3=0.1; beta4=0.05;

beta = [beta1,beta2,beta3,beta4];

A=100;

for r=1:A

sigma = 1;

epsilon = randn(n,1)*sigma;

y = beta1*X1+beta2*X2+beta3*X3+beta4*X4+epsilon;

b_est = X\y;

e = y - X*b_est;

sigma_est = (e'*e)/n;

V=sigma_est*inv(X'*X);

L= (((2*pi)/n)^(-n/2))*(exp(-n/2))*(e'*e)^(-n/2);

logL = -(n/2)*log(2*pi) - (n/2)*log(sigma_est) - (1/(2*sigma_est))*(e'*e);

%%Restricted model

beta3=0;

y_r = beta1*X1+beta2*X2+beta4*X4+epsilon;

X_r=[X1,X2,X4];

b_est_r=X_r\y_r;

e_r = y_r - X_r*b_est_r;

sigma_est_r = (e_r'*e_r)/n;

L_r= (((2*pi)/n)^(-n/2))*(exp(-n/2))*(e_r'*e_r)^(-n/2);

logL_r= -(n/2)*log(2*pi) - (n/2)*log(sigma_est_r) - (1/(2*sigma_est_r))*(e_r'*e_r);

%LM %%% not sure whether this is the correct way of calculation

s1=(1/sigma_est_r)*X_r'*e_r;

s2=-n/(2*sigma_est_r)+(1/(2*(sigma_est_r^2)))*(e_r'*e_r);

s=[s1; s2];

V_r=sigma_est_r*inv(X_r'*X_r);

LM(r)=s1'*V_r*s1;

end

Am I doing the LM calculation correctly? Thank you.