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
#include<oxstd.h>
#include<oxfloat.h>
#include<oxprob.h>
#include<oxdraw.h>
#include<quadpack.h>
#import<maximize>
//these variables have to be global variables as they are used in the functions
decl mMO, mLO;
decl dthreshold;
decl dmax;
decl integral;
decl par1, par2, par3;
decl abserror;
//function used inside the integral
fMOProb(x)
{
x = probexp(x - dthreshold, par3) *
((1/par2) .* densn((x - par1) / par2));
return 1, x;
}
//Loglikelihood function is defined
fLogLik(const vP, const adFunc, const avScore,
const amHessian)
{
par1 = vP[0];
par2 = vP[1];
par3 = vP[2];
//limit orders smaller than last transaction
decl vLOlb = selectifr(mLO, mLO .<= dthreshold);
decl vLOhb = selectifr(mLO, mLO .> dthreshold);
//first part of the loglik function for limit orders smaller than last transaction
decl dFirst = rows(vLOlb) * (-0.5 * log(M_2PI) - 2 * log(vP[1]))
- 0.5 * sumsqrc((vLOlb - vP[0])./vP[1]);
//second part for limit orders bigger than last transaction
decl dSecond = rows(vLOhb) * (- 0.5 * log(M_2PI) - 2 * log(vP[1]))
- vP[2] * sumc(vLOhb - dthreshold)
- 0.5 * sumsqrc((vLOhb - vP[0])./vP[1]);
//third part for market orders with integral
QAGI(fMOProb, dthreshold, 1, &integral, &abserror);
decl dThird = rows(mMO) * log(integral);
adFunc[0] = dFirst + dSecond + dThird;
return 1, integral; //returns 1 if succeeded
}
main()
{
decl fic = fopen("F:/Diplom/0106/buysidemax/output.log", "la");
//datamatrix is imported
decl mPV = loadmat("F:/Diplom/0106/buy/p&v0106b.mat");// |
loadmat("F:/0405/buy/p&v0405bpre.mat");
//value for cutting the limit orders
decl d1 = 54;
//market and limit orders are separated
mMO = selectifr(mPV[][0], mPV[][0] .== 0);
mLO = selectifr(mPV[][0], mPV[][0] .> d1);
mLO = sortc(mLO);
//the best price of the buy side is determined
dmax = maxc(mLO);
//the vector of different thresholds for the grid search is determined
decl vthreshold = range(dmax-8, dmax-0.1, 0.1)';
decl vtheta = <60; 1; 1>, dfunc;
//starting values
decl i;
decl vresults = new matrix[rows(vthreshold)][1];
//vector for holding the function values
//loop for the grid search begins here
for(i = 0; i <= rows(vthreshold)-1; i++)
{
dthreshold = vthreshold[i];
print("\n-------------------------------------------------");
print("\nloop has arrived at i = ", i+1," of ", rows(vthreshold));
decl est = MaxBFGS(fLogLik, &vtheta, &dfunc, 0, TRUE);
vresults[i] = dfunc;
}
vthreshold = vresults ~ vthreshold;
//the 'best' threshold value is determined
dthreshold = minc(selectifr(vthreshold[][1], vthreshold[][0] .==
maxc(vthreshold[][0])));
print("\nThreshold value corresponding to the highest Likelihood: ", dthreshold);
//moments of limit orders are printed
print("\n", "\n The moments of the limit orders are: \n", "%r",
{"SampleSize", "Mean", "SD", "Skewness", "Kurtosis"}, moments(mLO));
//maximization is done with respect to the 'best' threshold
decl ir = MaxBFGS(fLogLik, &vtheta, &dfunc, 0, TRUE);
//The numerical hessian matrix is computed
decl mHessian;
Num2Derivative(fLogLik, vtheta, &mHessian);
//The numerical gradient is computed
decl vgradient;
Num1Derivative(fLogLik, vtheta, &vgradient);
//covariance matrix computed over the outer product of gradients
decl mCovOG = invertgen(vgradient * vgradient');
decl mSEOG = sqrt(mCovOG[0][0]~mCovOG[1][1]~mCovOG[2][2]);
decl mTOG = (vtheta .* (1 ./ mSEOG)');
decl mPOG = (tailt(sqrt(mTOG .^ 2), 2));
//covariance-matrix, standarderrors, t-values and p-values are computed
decl mCov = 1/(rows(mLO)) .* (-invertgen(mHessian, 0));
decl mSE = sqrt(mCov[0][0]~mCov[1][1]~mCov[2][2]);
decl mT = (vtheta .* (1 ./mSE)');
decl mP = (tailt(sqrt(mT .^ 2), 2));
//results are printed into the output window
print("\n Function value is ", dfunc, "\n at parameter values: ",
"%c", {"mu", "sigma", "lambda"}, "%r",
{"coeff", "StandErr", "t-Value", "p-Value", "StandardErr(OG)", "t-
value(OG)", "p-Value"},
vtheta'|mSE|mT'|mP'|mSEOG|mTOG'|mPOG' , "\n",
MaxConvergenceMsg(ir));
println("\nThe MOs probability is: ", integral, "\nwith absolute error: ", abserror,
"\n\nEstimated Number of MOs: ", integral * rows(mPV),
"\n\nTrue Number of MOs: ", rows(mMO));
print("\n\nThe covariance-matrix is: ", mCov);
print("\n\nThe covariance-matrix from outer products of gradients is: ",
mCovOG);
//sequences to draw the density are determined
decl u = range(minc(mLO)-2, dthreshold, 0.1)';
decl y = range(dthreshold, dmax+4, 0.1)';
//the densities for drawing are declared
decl dens1 = invert(vtheta[1]) .* densn((u-vtheta[0])./vtheta[1]);
decl dens2 = (1 - probexp(y - dthreshold, vtheta[2])) .*
(invert(vtheta[1]) .* densn((y - vtheta[0])./vtheta[1]));
decl dens3 = invert(vtheta[1]) .* densn((y-vtheta[0])./vtheta[1]);
SetDraw(SET_COLORMODEL, 4);
DrawXMatrix(0, vthreshold[][0]', "Loglikelihood", vthreshold[][1]', "Thresholds", 0,
2);
SaveDrawWindow("F:/Diplom/0106/buysidemax/gridsearch0106b.ps");
CloseDrawWindow(); //window with the grid search graphic is closed
DrawAdjust(ADJ_AREA_X, 0, minc(mLO)-2, dmax+4);
DrawAdjust(ADJ_AREAMATRIX, 2, 1, 1);
DrawXMatrix(1, log(selectifr(mPV[][1], mPV[][0] .!= 0 .&& mPV[][0] .> d1))',
"log(Volumes)", selectifr(mPV[][0], mPV[][0] .!= 0 .&& mPV[][0].>
d1)',"Prices(LO)", 1, 17);
DrawXMatrix(0, dens1', "Density", u', "Prices", 0, 2);
DrawXMatrix(0, dens2', "", y', "", 0, 2);
DrawXMatrix(0, dens3', "", y', "", 0, 21);
DrawXMatrix(0, constant(0, mLO)', "", mLO', "", 1, 3);
DrawDensity(0, mLO', "", 1, 0, 0, 0, 0, 0, 17);
DrawTitle(0, "Density of Prices (estimated by MSL)");
DrawTitle(1, "Scatterplot of log(Volumes) and Prices with all considered Limit
Orders");
SaveDrawWindow("F:/Diplom/0106/buysidemax/density0106b.ps");
}
