[161] | 1 | #include <estim/arx.h> |
---|
| 2 | #include <stat/libEF.h> |
---|
| 3 | #include <stat/loggers.h> |
---|
| 4 | using namespace itpp; |
---|
| 5 | |
---|
| 6 | //These lines are needed for use of cout and endl |
---|
| 7 | using std::cout; |
---|
| 8 | using std::endl; |
---|
| 9 | |
---|
| 10 | int main() { |
---|
| 11 | // Setup model |
---|
| 12 | RV y("1","{y }","1","0"); |
---|
| 13 | RV u1("2","{u1 }","1","0"); |
---|
| 14 | RV u2("3","{u2 }","1","0"); |
---|
| 15 | |
---|
| 16 | // Full system |
---|
| 17 | vec thg("0.7 1 1"); //Simulated system |
---|
| 18 | //y=a y_t-1 + u1 + u2 |
---|
| 19 | double sqr=0.1; |
---|
| 20 | int ord = 1; |
---|
| 21 | |
---|
| 22 | // Estimated systems ARX(2) |
---|
| 23 | RV thri("4","{thr }",vec_1(2+1),"0"); |
---|
| 24 | RV thrg("5","{thr }",vec_1(3+1),"0"); |
---|
| 25 | // Setup values |
---|
| 26 | |
---|
| 27 | //ARX constructor |
---|
| 28 | mat V0 = 0.001*eye(thri.count()); V0(0,0)*= 10; // |
---|
| 29 | mat V0g = 0.001*eye(thrg.count()); V0g(0,0)*= 10; // |
---|
| 30 | double nu0 = ord+1; |
---|
| 31 | double frg = 0.94; |
---|
| 32 | |
---|
| 33 | ARX P1(thri, V0, nu0, frg); |
---|
| 34 | ARX P2(thri, V0, nu0, frg); |
---|
| 35 | ARX PG(thrg, V0g, nu0, frg); |
---|
| 36 | //Test estimation |
---|
| 37 | int ndat = 1000; |
---|
| 38 | int t; |
---|
| 39 | |
---|
| 40 | // Logging |
---|
| 41 | dirfilelog L("exp/merg",1000); |
---|
| 42 | |
---|
| 43 | int Eth1_log = L.add(thri,"P1"); |
---|
| 44 | int Eth2_log = L.add(thri,"P2"); |
---|
| 45 | int Data_log = L.add(RV("10","{Y U1 U2 }","3","0"), ""); |
---|
| 46 | int LL_log = L.add(RV("11","{1 2 G }","3","0"), "LL"); |
---|
| 47 | |
---|
| 48 | vec Yt(ndat); |
---|
| 49 | |
---|
| 50 | Yt.set_subvector(0,randn(ord)); //initial values |
---|
| 51 | vec rgrg(thg.length()); |
---|
| 52 | vec rgri(thri.count()); |
---|
| 53 | |
---|
| 54 | for (t=2; t<ndat; t++) { |
---|
| 55 | // True system |
---|
| 56 | rgrg(0)=Yt(t-1); |
---|
| 57 | rgrg(1)= sin(t/(40)*pi); |
---|
| 58 | rgrg(2)= cos(t/(40)*pi); |
---|
| 59 | |
---|
| 60 | Yt(t) = thg*rgrg + sqr * NorRNG(); |
---|
| 61 | |
---|
| 62 | // 1st |
---|
| 63 | rgri(0)=Yt(t); |
---|
| 64 | rgri(1)=rgrg(1); |
---|
| 65 | P1.bayes(rgri); |
---|
| 66 | // 2nd |
---|
| 67 | rgri(1)=rgrg(2); |
---|
| 68 | P2.bayes(rgri); |
---|
| 69 | |
---|
| 70 | //Global |
---|
| 71 | PG.bayes(concat(Yt(t),rgrg)); |
---|
| 72 | L.logit(Eth1_log, P1._epdf().mean()); |
---|
| 73 | L.logit(Eth2_log, P2._epdf().mean()); |
---|
| 74 | L.logit(Data_log, vec_3(Yt(t), rgrg(1), rgrg(2))); |
---|
| 75 | L.logit(LL_log, vec_3(P1._ll(), P2._ll(), PG._ll())); |
---|
| 76 | L.step(false); |
---|
| 77 | } |
---|
| 78 | L.step(true); |
---|
| 79 | |
---|
| 80 | } |
---|