root/mpdm/merg_pred.cpp @ 161

#include <estim/arx.h>
#include <stat/libEF.h>
#include <stat/loggers.h>
using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main() {
        // Setup model
        RV y("1","{y }","1","0");
        RV u1("2","{u1 }","1","0");
        RV u2("3","{u2 }","1","0");
        
        // Full system
        vec thg("0.7 1 1"); //Simulated system 
                                                //y=a y_t-1 + u1 + u2
        double sqr=0.1;
        int ord = 1;
        
        // Estimated systems ARX(2)
        RV thri("4","{thr }",vec_1(2+1),"0");
        RV thrg("5","{thr }",vec_1(3+1),"0");
        // Setup values
        
        //ARX constructor
        mat V0 = 0.001*eye(thri.count()); V0(0,0)*= 10; //
        mat V0g = 0.001*eye(thrg.count()); V0g(0,0)*= 10; //
        double nu0 = ord+1;
        double frg = 0.94;
        
        ARX P1(thri, V0, nu0, frg);
        ARX P2(thri, V0, nu0, frg);
        ARX PG(thrg, V0g, nu0, frg);
        //Test estimation
        int ndat = 1000;
        int t;
        
        // Logging
        dirfilelog L("exp/merg",1000);

        int Eth1_log = L.add(thri,"P1");
        int Eth2_log = L.add(thri,"P2");
        int Data_log = L.add(RV("10","{Y U1 U2 }","3","0"), "");
        int LL_log   = L.add(RV("11","{1 2 G }","3","0"), "LL");
        
        vec Yt(ndat);
        
        Yt.set_subvector(0,randn(ord)); //initial values
        vec rgrg(thg.length());
        vec rgri(thri.count());
        
        for (t=2; t<ndat; t++) {
                // True system
                rgrg(0)=Yt(t-1); 
                rgrg(1)= sin(t/(40)*pi);
                rgrg(2)= cos(t/(40)*pi);
                
                Yt(t) = thg*rgrg + sqr * NorRNG();
                
                // 1st
                rgri(0)=Yt(t);
                rgri(1)=rgrg(1);
                P1.bayes(rgri);
                // 2nd
                rgri(1)=rgrg(2);
                P2.bayes(rgri);
                
                //Global
                PG.bayes(concat(Yt(t),rgrg));
                L.logit(Eth1_log, P1._epdf().mean());
                L.logit(Eth2_log, P2._epdf().mean());
                L.logit(Data_log, vec_3(Yt(t), rgrg(1), rgrg(2)));
                L.logit(LL_log, vec_3(P1._ll(), P2._ll(), PG._ll()));
                L.step(false);
        }
        L.step(true);
        
}