MaxLikCode

This Section contains a number of codes that are used in the Maximum Likelihood Section.

MLse.m

This code simulates a linear model, estimates it by PLS and ML and calculates standard errors. For this code to run you need to have the following function accessible to MATLAB (i.e. in the same folder or on a pre-specified search path: OLSest, nll_lin, HessMp, gradp. You will also need the optimization toolbox. If you do not have access to the optimisation toolbox you can replace fminunc with fminsearch as the later is part of the main MATLAB software package.

% Code to
% a) simulate linear model, y_i = 0.2 + 0.6 x_1i -1.5 x_2i + err_i
%       x_1i and x_s1 come from N(0,1)
%       err_i ~ N(0,sd^2)
% b) estimate it by OLS
% c) estimate ot by ML
% d) estimate different ML standard errors
%
% This code requires the following Functions:
% OLSest, nll_lin, HessMp, gradp
clc
clear all

T   = 1000;                     % set sample size
b0  = [0.2; 0.6; -1.5];         % set parameter values
sd  = 0.5;                      % set error standard deviation
x   = [ones(T,1) randn(T,2)];   % simulate X matrix
err = randn(T,1)*sd;            % simulate error terms
y = x*b0 + err;                 % calculate y_i s

[b,bse,res,n,rss,r2] = OLSest(y,x,1);   % OLS estimation

datamat = [y x];                        % define data matrix for use in nll_lin
theta0 = [mean(y); 0; 0; std(y)];       % this sets the initial parameter vector
options = optimset; % sets optimisation options to default
[thetaopt] = fminsearch(@nll_lin,theta0,options,datamat,1);
 
H = HessMp(@nll_lin,thetaopt,datamat,1); % this returns the negative of the Hessian
g = gradp(@nll_lin,thetaopt,datamat,0);  % this returns a (T x size(theta0,1)) matrix of gradients
J = (g'*g);                             % calculates the OPG

se_H = sqrt(diag(inv(H)));
se_J = sqrt(diag(inv(J)));
se_SW = sqrt(diag(inv(H*inv(J)*H)));    % Sandwich variance covariance

disp('   Est     se(OLS)     se(H)     se(J)     se(SW)');
disp([thetaopt [bse;0] se_H se_J se_SW]);

nll_lin.m

This is the negative log likelihood function for a linear model. Save this as nll_lin.m.

function nll = nll_lin( theta, data , vec)
% input:    (i) theta, coef vector, last element is error sd
%           (ii), data matrix, col1: y cols2:end: explan. variables (incl constant)
%           (iii), 0 = if vector of loglikelihoods, and 1 if sum should be
%           returned
beta = theta(1:end-1);
sig  = abs(theta(end))+0.000001;    % this ensures a non-zero variance
y = data(:,1);
x = data(:,2:end);
res = y - x*beta;
nll = (-0.5)*(-log(2*pi)-log(sig^2)-(res.^2/(sig^2)));

if vec
    nll = sum(nll);
end

end

Gradient and Hessian code

These two functions are needed in order to calculate the Hessian and Gradient.

HessMp.m

Download the file from here. An example call to this function can be seen in the above MLse.m code. The first input is a handle to the Function that is used to calculate the Hessian (here nll_lin). All other inputs follow the inputs required for that function (here (theta, data , vec)). When used for calculating the variance-covariance matrix of parameter estimates we want to feed in the estimated parameters (thetaopt). The Hessian routine expects the function to return a scalar likelihood function and hence the input for vec is 1 (see definition of nll_lin.m).

H = HessMp(@nll_lin,thetaopt,datamat,1); % this returns the negative of the Hessian

gradp.m

Download the file from here. An example call to this function can be seen in the above MLse.m code. The first input is a handle to the Function that is used to calculate the Gradient (here nll_lin). All other inputs follow the inputs required for that function (here (theta, data , vec)). When used for calculating the variance-covariance matrix of parameter estimates we want to feed in the estimated parameters (thetaopt). The Gradient will calculate the gradient for either the scalar log-likelihood (if the vec input is set to 1) or the gradient at every observation (if the vec input is set to 0 - see definition of nll_lin.m). To calculate the OPG we need the latter.

g = gradp(@nll_lin,betaopt,datamat,0);  % this returns a (T x size(theta0,1)) matrix of gradients