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