%% Assess Model Specifications Using the Lagrange Multiplier Test
% Compare two model specifications for simulated education and income data. The unrestricted model has the following loglikelihood:
% l(\beta,\rho) = -n\log\Gamma(\rho) + \rho\sum^n_{k=1}\log\beta_k +
% (\rho-1)\sum^n_{k=1}\log(y_k)-\sum^n_{k=1}y_k\beta_k,
% where
% \beta_k=\frac{1}{\beta+x_k}.
% x_k is the number of grades that person k completed.
% y_k is the income (in thousands of USD) of person k.
% That is, the income of person k given the number of grades that person k completed is Gamma distributed with shape \rho and rate \beta_i. The restricted model sets \rho = 1, which implies that the income of person k given the number of grades person k completed is exponentially distributed with mean \beta + x_i.
% The restricted model is H_0: \rho = 1. In order to compare this model to the unrestricted model, you require:
% The gradient vector of the unrestricted model
% The maximum likelihood estimate (MLE) under the restricted model
% The parameter covariance estimator evaluated under the MLEs of the restricted model
% Load the data.

load Data_Income1
x = DataTable.EDU;
y = DataTable.INC;

%
% Estimate the restricted model parameters by maximizing l(\rho,\beta) with respect to \beta subject to the restriction \rho = 1. The gradient of l(\rho,\beta) is
% \frac{\partial l(\rho,\beta)}{\partial \beta} = \sum^T_{i=1}(y_i\beta_i^2-\rho\beta_i)
% \frac{\partial l(\rho,\beta)}{\partial \rho} = -T\Psi(\rho)+\sum^T_{i=1}(\log\beta_iy_i),
% where \Psi(\rho) is the digamma function.

rho0 = 1; % Restricted rho
dof = 1;  % Number of restrictions
dLBeta = @(beta) sum(y./((beta + x).^2) - rho0./(beta + x));...
    % Anonymous gradient function

[betaHat0,fVal,exitFlag] = fzero(dLBeta,0)

beta = [0:0.1:50];
plot(beta,arrayfun(dLBeta,beta))
hold on
plot([beta(1);beta(end)],zeros(2,1),'k:')
plot(betaHat0,fVal,'ro','MarkerSize',10)
xlabel('{\beta}')
ylabel('Loglikelihood Gradient')
title('{\bf Loglikelihood Gradient with Respect to \beta}')
hold off

% The gradient with respect to \beta (dLBeta) is decreasing, which suggests that there is a local maximum at its root. Therefore, betaHat0 is the MLE for the restricted model.  fVal indicates that the value of the gradient is very close to 0 at betaHat0. The exit flag (exitFlag) is 1, which indicates that fzero found a root of the gradient without a problem.
%
% Estimate the parameter covariance under the restricted model using the outer product of gradients (OPG).

rGradient = [-rho0./(betaHat0+x)+y.*(betaHat0+x).^(-2),...
      log(y./(betaHat0+x))-psi(rho0)];    % Gradient per unit
rScore = sum(rGradient)';                 % Score function
rEstParamCov = inv(rGradient'*rGradient); % Parameter covariance estimate

%
% Test the unrestricted model against the restricted model using the Lagrange multiplier test.

[h,pValue] = lmtest(rScore,rEstParamCov,dof)

% pValue is close to 0, which indicates that there is strong evidence to suggest that the unrestricted model fits the data better than the restricted model.
% Copyright 2015 The MathWorks, Inc.