p74 probit4.ox

时间序列分析中著名的OxMetrics软件包

#include <oxstd.h>
#import <database>
#import <maximize>

eprint
enum
{   Y_VAR, X_VAR
};

/*-------------------------- Probit : Database -----------------------------*/
class Probit : Database
{
    decl m_iModel;                  /* TRUE if model successfully estimated */
    decl m_fPrint;               /* TRUE if output desired after estimation */
    decl m_iResult;                             /* return code from MaxBFGS */
    decl m_iT1est, m_iT2est;                           /* estimation sample */
    decl m_cT;                                    /* number of observations */
    decl m_cY;             /* number of dependent variables (always 1 here) */
    decl m_cX;                             /* number of regressors (if any) */

    decl m_mY;                                /* dependent variable [cT][1] */
    decl m_mX;                          /* regressor data vector [cT][m_cX] */
    decl m_mPar;                                      /* parameters: [p][1] */
    decl m_dLoglik;                                       /* log-likelihood */

    Probit();                                                /* constructor */
    InitData();                               /* extract data from database */
    Estimate(const vStart);                            /* do the estimation */
    Output();                                          /* print the results */
    fProbit(const vP, const adFunc, const avScore, const amHessian);  /* LL */
};
/*------------------------ END Probit : Database ---------------------------*/

/*-------------------------- Probit : Database -----------------------------*/
Probit::Probit()
{
    this->Database();           // intialize base class
    m_iModel = 0;               // 1 when estimated
    m_fPrint = TRUE;
    m_iResult = -1;

    m_mPar = 0;

    print("Probit class example 4, object created on ",
        date(), ".\n\n");
}

Probit::fProbit(const vP, const adFunc, const avScore,
    const amHessian)
{
    decl prob = probn(m_mX * vP);
    decl tail = 1 - prob;

    adFunc[0] = double(
        meanc(m_mY .* log(prob) + (1-m_mY) .* log(tail)));

    if (avScore)
    {
        decl weight = (m_mY - prob) .* densn(m_mX * vP)
            ./ (prob .* tail);
        avScore[0] = meanc(weight .* m_mX)';
    }

return 1;                           // 1 indicates success
}

Probit::InitData()
{
    decl cp, vp, mh;

    m_iT1est = m_iT1sel;  m_iT2est = m_iT2sel;

    m_mY = GetGroupLag(Y_VAR, 0, 0);
    m_cY = columns(m_mY);

    if (m_cY != 1)
    {   eprint("Only univariate models allowed\n");
        return FALSE;
    }
    m_mX = GetGroup(X_VAR);
    m_cX = columns(m_mX);
    if (!m_cX)
    {   eprint("Need some regressors\n");
        return FALSE;
    }
    m_cT = m_iT2est - m_iT1est + 1;
    if (m_cT <= 2)
    {   eprint("Only ", m_cT, " observations. This is not enough.\n");
        return FALSE;
    }

    return TRUE;
}

Probit::Estimate(const vStart)
{
    if (!InitData())
    {   print("Failed to load data\n");
        m_mPar = zeros(1, m_cX);
        return;
    }
    if (sizer(vStart))
        m_mPar = vStart;
    if (m_fPrint)
        print("starting values: ", m_mPar);

    m_iResult = MaxBFGS(fProbit, &m_mPar, &m_dLoglik, 0, FALSE);
    m_dLoglik *= m_cT;

    if (m_fPrint)
        Output();
}

Probit::Output()
{
    print("\n", MaxConvergenceMsg(m_iResult),
        " using analytical derivatives",
        "\nFunction value = ", m_dLoglik, "; parameters:", m_mPar);

                  // if converged: compute standard errors
    if (m_iResult == MAX_CONV || m_iResult == MAX_WEAK_CONV)
    {
        decl mhess;
        if (Num2Derivative(fProbit, m_mPar, &mhess))
        {
            mhess = -invert(mhess) / m_cT;
            print("standard errors:", sqrt(diagonal(mhess)'));
        }
    }
}

main()
{
    decl probitobj;

                       // create an object of class Probit
    probitobj = new Probit();

                             // load data file into object
    probitobj->LoadIn7("data/finney.in7");
    probitobj->Info();           // print database summary
    probitobj->Deterministic(FALSE);    // create constant

                                    // Formulate the model
    probitobj->Select(Y_VAR, { "vaso",0,0 } );
    probitobj->Select(X_VAR, { "Constant",0,0,
        "Lrate",0,0, "Lvolume",0,0 } );
    probitobj->SetSelSample(-1, 1, -1, 1);  // full sample

    MaxControl(-1, 1, 1);          // print each iteration
    probitobj->Estimate(<-0.465; 0.842; 1.439>);//maximize

    delete probitobj;
}