Dear all,
I am trying to implement an analytical gradient for the likelihood function L given in the attached pdf. I give the gradients I have calculated there as well.
My MATLAB implementation of the negative log-likelihood function and gradients is given below and I use fminunc to solve it. I have run the optimization on two different datasets (attached) and get differences in the analytical and numeric gradients on the order of 1.8e-6 on one dataset (set ex = 2) but not on the other (set ex = 1). Also my parameter estimates do not match those obtained from running equivalent regressions in STATA with differences around 5% of the parameter value on average.
I am assuming I have some error in the analytical gradient implementation of the beta-components of the gradient vector. If anyone can help me out here and spot the error I would be very grateful. I will of course provide additional information as needed. The function and script needed to reproduce the behavior are below.
Thanks very much,
Peter
Here is the implementation of the log likelihood function and gradients and the data files are attached.
function [negLogLik grad] = truncReg_cost(Z,theta,params,c)
sigma = sqrt(params(1));
beta = params(2:end);
nObs = size(Z,1);
c = ones(nObs,1)*c;
err = theta-Z*beta;
sse = err'*err;
cdfArg = (c-Z*beta)/sigma;
quotDens = normpdf(cdfArg,0,1)./(ones(nObs,1)-normcdf(cdfArg,0,1));
negLogLik = nObs*log(sigma*sqrt(2*pi))+0.5*(1/sigma^2)*sse+...
sum(log(ones(nObs,1)-normcdf(cdfArg,0,1)));
partialSigma = 0.5*((nObs/sigma^2)-(1/sigma^4)*sse+...
(1/sigma^3)*(cdfArg*sigma)'*quotDens);
partialBeta = -(1/sigma^2)*Z'*err+...
(1/sigma)*Z'*quotDens;
grad = [partialSigma;partialBeta];
end
To reproduce the issue you can run this script
ex = 1;
if ex == 1
data = csvread('trExample.csv');
else
data = csvread('trExample2.csv');
end
if ex == 1
c = 0;
truncDir = 'l';
indVar = [ones(size(data,1),1) data(:,3:end)];
depVar = data(:,2);
else
c = 40;
truncDir = 'l';
indVar = [ones(size(data,1),1) data(:,2:end)];
depVar = data(:,1);
end
if truncDir == 'l'
indVar = indVar(depVar>c,:);
depVar = depVar(depVar>c,:);
else
indVar = indVar(depVar<c,:);
depVar = depVar(depVar<c,:);
end
betaInit = pinv(indVar'*indVar)*indVar'*depVar;
sigmaInit = var(depVar);
params = [sigmaInit; betaInit];
costfunc = @(params) truncReg_cost(indVar,depVar,params,c);
options = optimset('GradObj','on','Algorithm','interior-point', ...
'Display','on','DerivativeCheck','on');
[outParams,LL,ef,out,grad] = fminunc(costfunc,params,options);