ConScore
==============================================

Purpose
----------------

Compute local score statistic and its probability for hypotheses involving parameters under constraints

Format
----------------
.. function:: { SL, SLprob } = ConScore(H, G, grad, a, b, c, d, bounds, psi)

    :param H: Hessian of loglikelihood with respect to parameters.
    :type H: KxK matrix

    :param G: cross-product matrix of the first derivatives by observation. If not available set to *H*.
    :type G: KxK matrix

    :param grad: gradient of loglikelihood with respect to parameters.
    :type grad: Kx1 vector

    :param a: linear equality constraint coefficients.
    :type a: MxK matrix

    :param b: linear equality constraint constants.

        These arguments specify the linear equality
        constraints of the following type:

        .. math::  a * X = b

        where *X* is the Kx1 parameter vector.
    :type b: Mx1 vector

    :param c: linear inequality constraint coefficients.
    :type c: MxK matrix

    :param d: linear inequality constraint constants.

        These arguments specify the linear inequality
        constraints of the following type:
        .. math::  c * X >= d

        where *X* is the Kx1 parameter vector.
    :type d: Mx1 vector

    :param bounds: bounds on parameters. The first column
        contains the lower bounds, and the second column the
        upper bounds.
    :type bounds: Kx2 matrix

    :param psi: indices of the set of parameters in the hypothesis
    :type psi: matrix

    :return SL: local score statistic of hypothesis.

    :rtype SL: scalar

    :return SLprob: probability of *SL*.

    :rtype SLprob: scalar

Examples
----------------
This example is from Silvapulle and Sen, Constrained Statistical Inference, page 181-3. It computes the
local score statistic and probability for an ARCH
model. It tests the null hypothesis of no arch effects against the alternative of arch effects
subject to their being constrained to be positive.
The Hessian, *H*, cross-product matrix, *G*, and the
gradient vector, *grad*, are generated by an
estimation using :func:`sqpSolvemt` where the model is
an ARCH model with the arch parameters constrained to be zero.

::

    // dataMatrix
    struct DS d0;
    d0 = reshape(dsCreate,2,1);

    load z0[] = aoi.asc;
    z = packr(lagn(251*ln(trimr(z0, 1, 0)./trimr(z0, 0, 1)), 0|1|2|3|4));
    d0[1].dataMatrix = z[., 1];
    d0[2].dataMatrix = z[., 2:5];

    // Control structure
    struct sqpsolvemtControl c0;
    c0 = sqpSolveMTcontrolCreate;

    /*
    ** Constraints setting arch parameter equal
    ** to zero for H(theta) = 0
    */
    c0.A = zeros(3, 6) ~ eye(3);
    c0.B = zeros(3, 1);

    /*
    ** Causes cross-product of
    ** Jacobian to be computed which
    ** is needed for ConScore
    */
    c0.covType = 2;

    struct PV p0;
    p0 = pvPack(pvCreate, .08999, "constant");
    p0 = pvPack(p0, .25167|-.12599|.09164|.07517,
        "phi");
    p0 = pvPack(p0,3.22713, "omega");
    p0 = pvPack(p0, 0|0|0, "arch");


    struct sqpsolvemtOut out0;
    out0 = sqpsolvemt(&lpr, p0, d0, c0);

    // Set up constraints for H(theta) >= 0
    bounds = { -1e256 1e256,
               -1e256 1e256,
               -1e256 1e256,
               -1e256 1e256,
               -1e256 1e256,
               -1e256 1e256,
                    0 1e256,
                    0 1e256,
                    0 1e256 };
    H = out0.hessian;
    G = out0.xproduct;

    // minus because -logl in log-likelihood
    grad = -out0.gradient;

    psi = { 7, 8, 9 };

    { SL, SLprob } = ConScore(H,G,grad,0,0,0,0,bounds,psi);

will assign the variables *SL* and *SLprob* as follows:

::

    SL = 3.8605086  
    SLprob = 0.10410000

Remarks
-------

:func:`ConScore` computes the local score statistic for the hypothesis :math:`H(Θ) = 0`
vs. :math:`H(Θ) ≥ 0`, where :math:`Θ` is the vector of estimated parameters, and :math:`H()` is
a constraint function of the parameters.

First, the model with :math:`H(Θ) = 0` is estimated, and the Hessian and
optionally the cross-product of the derivatives is computed. Also, the
gradient vector is computed.

Next, the constraint arguments are set to :math:`H(Θ) ≥ 0`.