Dear Simon,
Here is the function "ScoreContributions" due to Michael Creel, that
does the job.
I trust he won't mind me passing it on here bundled with a companion
"HessianContributions" function, which does the second derivatives,
should you need them.
Note that "func" should match the format as in Num1Derivative except
that the second argument should pass the address of a column vector,
not a double.
Best
James
----------------------------------------------------------------
James Davidson
Professor of Econometrics
School of Business and Economics, University of Exeter
TSM Email: [log in to unmask]
TSM Web: http://www.timeseriesmodelling.com
----------------------------------------------------------------
/*--------------------------------ScoreContributions
------------------------------*/
/*
This is a modification of Num1Derivative that allows the
function to
be vector valued, as with Gauss's gradp. The function returns n
x 1,
the parameter is p x 1, and the derivative is n x p. The matrix
of
derivatives is passed as the third argument, just as with
Num1Derivative.
Michael Creel, [log in to unmask]
21 Aug. 1998
*/
const decl SQRT_EPS =1E-8; /* appr. square root of machine
precision */
const decl DIFF_EPS1=5E-6; /* Rice's formula:
log(DIFF_EPS)=log(MACH_EPS)/3 */
dFiniteDiff1(const x)
{
return max( (fabs(x) + SQRT_EPS) * SQRT_EPS, DIFF_EPS1);
}
//**********************************************************/
ScoreContributions(const func, vP, const avScore)
{
decl i, cp = rows(vP), left, right, fknowf = FALSE, p, h, f, fm, fp,
v;
for (i = 0; i < cp; i++) /* get 1st derivative by central
difference */
{
p = double(vP[i][0]);
h = dFiniteDiff1(p);
vP[i][0] = p + h;
right = func(vP, &fp, 0, 0);
if(i==0)
v = new matrix[rows(fp)][cp];
vP[i][0] = p - h;
left = func(vP, &fm, 0, 0);
vP[i][0] = p; /* restore original
parameter */
if (left && right)
v[][i] = (fp - fm) / (2 * h); /* take central
difference */
else
return FALSE;
}
avScore[0] = v;
return TRUE;
}
//*********************************************************
HessianContributions(const func, vP, const avHess)
{
decl cp = rows(vP), leftleft, leftright, rightleft, rightright,
p, q, h, fpp, fpq, fqp, fqq;
decl v = <>;
for (decl i = 0; i < cp; i++)
for (decl j = i; j < cp; j++)
{
p = double(vP[i][0]);
q = double(vP[j][0]);
h = dFiniteDiff1(p);
vP[i][0] += h;
vP[j][0] += h;
rightright = func(vP, &fpp, 0, 0);
vP[i][0] -= 2*h;
leftright = func(vP, &fqp, 0, 0);
vP[j][0] -= 2*h;
leftleft = func(vP, &fqq, 0, 0);
vP[i][0] += 2*h;
rightleft = func(vP, &fpq, 0, 0);
vP[i][0] = p; /* restore original
parameter */
vP[j][0] = q;
if ((leftleft && leftright) && (rightleft && rightright))
v ~= (fpp - fpq - fqp + fqq) / (4 * h^2); /* take
central difference */
else return FALSE;
}
avHess[0] = v;
return TRUE;
}
/*----------------------End of ScoreContributions
------------------------------*/
> -----Original Message-----
> From: The ox-users list is aimed at all Ox users
> [mailto:[log in to unmask]] On Behalf Of Simon
> Sent: 31 July 2008 10:19
> To: [log in to unmask]
> Subject: Covariance Estimation with outer product of gradients
>
> Dear Ox-Users,
>
>
> I use Ox for now more than a year successfully, but with this
> problem I do not know what to do. I am pretty sure that
> someone of you professionals can give me an advice how to proceed.
>
> The following:
>
> I have programed an Ox-Code to estimate a very specific
> Maximum Likelihood function with an integral in it which is
> approximated via Quad-package. Everything works quite well. I
> get my parameter estimates and everything what I want to. The
> code is allright, but
>
>
> I want to estimate the covariance of my estimators by the
> outer product of gradients. The gradient given out via
> Num1Derivative() is the value of the numerical first
> derivative at my estimated values for the whole likelihood
> function right? So it is nearly zero resulting in extremely
> high standard errors.
>
> What I need is the gradient for every single observation in
> the loglikelihood function.
> Such that I can build for every observation an outer product
> of gradients and sum the resulting matrices. I cannot see any
> possibility how to proceed. Can anybody help me with good ideas?
>
> I provided the code below and hope that it works to attach
> the data as a .mat-file.
>
> Thank you all!
>
>
> Simon
|