#include <oxstd.h>

#import <pcnaive>
#import <pcfiml>
#import <simula>


/*---------------- SimSystem : Simulation ---------------*/
class SimSystem : Simulation
{
    decl m_dgp;                            // DGP container
    decl m_sys;                         // system container
    decl m_mCoef;                           // coefficients
    decl m_mTest;                        // test statistics
    SimSystem(const mxT);                    // constructor
    ~SimSystem();                             // destructor
};
SimSystem::SimSystem(const mxT)
{
    m_dgp = new PcNaiveDgp(3,3);          // create the DGP

    m_sys = new PcFiml();              // create the system
    m_sys.Create(1, 1, 1, mxT, 1);   //create the database
    m_sys.Append(zeros(mxT, 6),         // initialize data
        {"YA", "YB", "YC", "ZA", "ZB", "ZC"}, 0);
    m_sys.Deterministic(FALSE);   // create constant,trend

    m_sys.SetPrint(FALSE);       // no printing of results
}
SimSystem::~SimSystem()
{
    delete m_dgp;                         // delete the dgp
    delete m_sys;                         // and the system
}
/*----------- END SimSystem : Simulation ----------------*/

/*-------------------- SimArTest1 -----------------------*/
class SimArTest1 : SimSystem
{
    decl m_fUseZ;            // TRUE if Z used in dgp/model
    decl m_fUseYlag;       // TRUE if lagged Y in dgp/model
    decl m_cTdiscard;      // no of observations to discard

    SimArTest1(const cRep, const mDgpY_1, const mDgpZ,
        const mDgpEvar, const mDgpE_ar);     // constructor
    ~SimArTest1();                            // destructor
    Generate(const iRep, const cT,  // generate replication
        const mxT);
    GetTestStatistics();          // return test statistics
    GetPvalues();               // return p-values of tests
};
SimArTest1::SimArTest1(const cRep, const mDgpY_1, const mDgpZ,
    const mDgpEvar, const mDgpE_ar)
{
    decl mt, mxt, z31 = zeros(3,1), z33 = zeros(3,3);

    mt = <50,100>;          // sample sizes for Monte Carlo
    m_cTdiscard = 20;                 // use 20 to start up
    mxt = m_cTdiscard + max(mt);          // maximum needed

    this.SimSystem(mxt);          // call base constructor

    m_fUseZ = any(mDgpZ .!= 0); // check if Z in Y equation
    m_fUseYlag = any(mDgpY_1 .!= 0); //or lagged Y in Y eqn

                                              // set up DGP
    m_dgp.SetYParameter(z33, mDgpY_1, mDgpZ, z31);
    m_dgp.SetUParameter(mDgpE_ar, z33);
    m_dgp.SetZParameter(mDgpEvar, z31, z31);
    m_dgp.SetDistribution(U_DGP, MVNORMAL, z31, mDgpEvar);

                                     // formulate the model
    m_sys.Select(Y_VAR,{ "YA",0,0, "YB",0,0, "YC",0,0} );
    m_sys.Select(X_VAR,{ "Constant",0,0 } );

    if (m_fUseZ)           // model has Z vars, as does DGP
        m_sys.Select(X_VAR,{"ZA",0,0,"ZB",0,0,"ZC",0,0});
    if (m_fUseYlag)      // model has lagged Y, as does DGP
        m_sys.Select(Y_VAR,{"YA",1,1,"YB",1,1,"YC",1,1});

                               // set simulation parameters
    this.Simulation(mt, mxt, cRep, TRUE, -1,
        <0.2,0.1,0.05,0.01>, 0);

    this.SetTestNames({"normal", "F AR 1", "Chi AR 1"});

//  m_sys.Output(TRUE, FALSE);
    m_dgp.Print();
//  m_dgp.Asymp();

    print(m_cTdiscard, " observations discarded\n");
}
SimArTest1::~SimArTest1()
{
    ~SimSystem();                   // call base destructor
}
SimArTest1::Generate(const iRep, const cT, const mxT)
{
    decl i, y, test;

    if (iRep && imod(iRep,10000) == 0)
        eprint("Repetition ",iRep," sample size ",cT,"\n");
    if (iRep == 0)
        m_dgp.DiscardZ();                    // need new Z

    y = m_dgp.Generate(mxT);          // generate the data

    m_sys.Renew(y[][0], "YA", 0);     // store in database
    m_sys.Renew(y[][1], "YB", 0);
    m_sys.Renew(y[][2], "YC", 0);

    if (iRep == 0)
    {   if (m_fUseZ)                  // yes: also update Z
        {   decl z = m_dgp.GetZ();
            m_sys.Renew(z[][0], "ZA", 0);
            m_sys.Renew(z[][1], "ZB", 0);
            m_sys.Renew(z[][2], "ZC", 0);
        }
        m_sys.SetSample(m_cTdiscard + 1, 1,
            m_cTdiscard + cT, 1); // set the correct sample
    }
    m_sys.Estimate();                // estimate the model
                  // generate and store the test statistics
    m_mTest = m_sys.NormalityTest() ~ m_sys.ArTest(1,1);

return TRUE;
}
SimArTest1::GetTestStatistics()
{
    return m_mTest[0][]';
}
SimArTest1::GetPvalues()
{
    return m_mTest[1][]';
}
/*------------------- END SimArTest1 --------------------*/

main()
{
    decl exp, gamma0, gamma1, gamma2, z33, o33;

    gamma0 = ones(3,3) / 10 + unit(3) / 5;
    gamma1 = <0.5,0.1,0;0.1,0.5,0;0,0,0.3>;
    gamma2 = <0.5,0.3,0;0.3,0.7,0;0,0,0.7>;
    z33 = zeros(3,3);
    o33 = ones(3,3);

    eprint("(a): white noise (size)\n");
    print("\n(a): white noise (size) =============\n");
	decl time = timer();
    exp = new SimArTest1(10000, z33, z33, gamma1, z33);
    exp.Simulate();
    delete exp;
	println("time: ", timespan(time));
/*
    eprint("(b): no dynamics (size)\n");
    print("\n(b): no dynamics (size) =============\n");
    exp = new SimArTest1(1000, z33, o33, gamma1, z33);
    exp.Simulate();
    delete exp;

    eprint("(c): VAR(1) (size)\n");
    print("\n(c): VAR(1) (size) =============\n");
    exp = new SimArTest1(1000, gamma0, z33, gamma1, z33);
    exp.Simulate();
    delete exp;

    eprint("(e): VARX(1) (size)\n");
    print("\n(e): VARX(1) (size) =============\n");
    exp = new SimArTest1(1000, gamma0, o33, gamma1, z33);
    exp.Simulate();
    delete exp;

    eprint("(g): autocorrelation\n");
    print("\n(g): autocorrelation =============\n");
    exp = new SimArTest1(1000, gamma0, z33, gamma1, gamma1);
    exp.Simulate();
    delete exp;

    eprint("(h): autocorrelation\n");
    print("\n(h): autocorrelation =============\n");
    exp = new SimArTest1(1000, gamma0, z33, gamma1, gamma2);
    exp.Simulate();
    delete exp;
*/
}
