#include <oxstd.h>
// Original: #include <oxdraw.h>
#include <packages/gnudraw/gnudraw.h>

#include <simula.h>

#pragma  link("simula.oxo")

/*-------------------------- test statistics -------------------------------*/
/*------------------------- ENDtest statistics -----------------------------*/

/*------------------------- SimAR1 : Simulation ----------------------------*/
class SimAR1 : Simulation         // inherit from simulation class
{
    decl m_mTest;                              // test statistic
    decl m_mCoef;                              // coefficients
    decl m_cTdiscard;                          // no of obs to discard
    decl m_dAr;                                // AR parameter
    decl m_cDf;                                // degrees of freedom
    decl m_mAr1Val;                            // ar-values of each replication

    SimAR1(const dAr, const cT);               // constructor
    ~SimAR1();                                 // destructor
    Generate(const iRep, const cT, const mxT); // generate replication
    GetPvalues();                              // return p-values of tests
    GetCoefficients();                         // return coefficients
    GetTestStatistics();                       // return test statistics
    Plot(const iRep);                          // plot function
};

SimAR1::SimAR1(const dAr, const cT)
{
    decl crep = 2000;
    m_cTdiscard = 20;
    m_dAr = dAr;
    m_mAr1Val = zeros(1, crep);

    Simulation(matrix(cT), cT, crep, TRUE, -1, <0.2,0.1,0.05,0.01>,
        matrix(m_dAr));

    SetCoefNames({"ar-param"});
    SetTestNames({"t-value"});
}
SimAR1::~SimAR1()
{
}
SimAR1::Plot(const iRep)
{
    if (iRep >= 100 && imod(iRep + 1, 100) == 0)
    {
        DrawDensity(0, m_mAr1Val[0][:iRep - 1], "", FALSE, TRUE, TRUE);
        DrawQQ(1, m_mAr1Val[0][:iRep - 1] - m_dAr, "", QQ_N, m_cDf, 0);

//      DrawDensity(2, m_mTestVal[0][:iRep - 1], "", FALSE, TRUE, TRUE);
//      DrawQQ(3, m_mTestVal[0][:iRep - 1], "", QQ_T, m_cDf, 0);

        ShowDrawWindow();
    }
}
SimAR1::Generate(const iRep, const cT, const mxT)
{
    decl eps, y, ylag, xtxinv, ar1, res, sigmasqr;

    eps = rann(cT + m_cTdiscard, 1);                 /* generate error term */
    y = cumulate(eps, matrix(m_dAr));
    ylag = lag0(y, 1);

    y = y[m_cTdiscard:][];                 /* discard intitial observations */
    ylag = ylag[m_cTdiscard:][];

    xtxinv = invert(ylag'ylag);
    ar1 = (xtxinv * (ylag'y))[0][0];
    res = y - ylag * ar1;
    m_cDf = cT - 1;
    sigmasqr = double(res'res / m_cDf);
//print(ar1, " ", sqrt(sigmasqr * xtxinv[0][0]), " ", sigmasqr, "\n");
    Plot(iRep);

    m_mTest = (ar1 - m_dAr) / sqrt(sigmasqr * xtxinv);
    m_mCoef = ar1;
    m_mAr1Val[0][iRep] = ar1;

return 1;                                 /* 1 indicates success, 0 failure */
}
SimAR1::GetPvalues()
{
    return matrix( 2 * tailt( fabs(m_mTest[0][]), m_cDf ) );
}
SimAR1::GetTestStatistics()
{
    return matrix(m_mTest);
}
SimAR1::GetCoefficients()
{
    return matrix(m_mCoef);
}
/*----------------------- END SimAR1 : Simulation --------------------------*/

main()
{
    decl experiment = new SimAR1(0.5, 50);                    // create object
    SetDrawWindow("AR 1, 0.5, T=50");
    experiment->Simulate();                                  // do simulations
    delete experiment;                                        // remove object

    experiment = new SimAR1(1.0, 50);                         // create object
    SetDrawWindow("AR 1, 1.0, T=50");
    experiment->Simulate();                                  // do simulations
    delete experiment;                                        // remove object
}