[161] | 1 | #include <estim/arx.h> |
---|
[198] | 2 | #include <estim/merger.h> |
---|
[384] | 3 | #include <stat/exp_family.h> |
---|
[161] | 4 | #include <stat/loggers.h> |
---|
[176] | 5 | //#include <stat/merger.h> |
---|
[254] | 6 | using namespace bdm; |
---|
[161] | 7 | |
---|
| 8 | //These lines are needed for use of cout and endl |
---|
| 9 | using std::cout; |
---|
| 10 | using std::endl; |
---|
| 11 | |
---|
| 12 | int main() { |
---|
| 13 | // Setup model |
---|
[162] | 14 | RV y ( "{y }" ); |
---|
| 15 | RV u1 ( "{u1 }" ); |
---|
| 16 | RV u2 ( "{u2 }" ); |
---|
[213] | 17 | RV ym=y; ym.t ( -1 ); |
---|
| 18 | RV yy = y; yy.add ( ym ); |
---|
| 19 | RV uu=u1; uu.add ( u2 ); |
---|
[162] | 20 | |
---|
[161] | 21 | // Full system |
---|
[198] | 22 | vec thg ( "0.7 1 1 0" ); //Simulated system - zero for constant term |
---|
[162] | 23 | //y=a y_t-1 + u1 + u2 |
---|
[161] | 24 | double sqr=0.1; |
---|
| 25 | int ord = 1; |
---|
[162] | 26 | |
---|
[161] | 27 | // Estimated systems ARX(2) |
---|
[198] | 28 | RV thri ( "{thr_i }",vec_1 ( 3+1 ) ); |
---|
| 29 | RV thrg ( "{thr_g }",vec_1 ( 4+1 ) ); |
---|
[161] | 30 | // Setup values |
---|
[162] | 31 | |
---|
[161] | 32 | //ARX constructor |
---|
[270] | 33 | mat V0 = 0.001*eye ( thri._dsize() ); V0 ( 0,0 ) *= 10; // |
---|
| 34 | mat V0g = 0.001*eye ( thrg._dsize() ); V0g ( 0,0 ) *= 10; // |
---|
[198] | 35 | double nu0 = ord+6.0; |
---|
[204] | 36 | double frg = 0.95; |
---|
[162] | 37 | |
---|
[270] | 38 | ARX P1; P1.set_statistics( 2, V0, nu0); P1.set_parameters(frg ); |
---|
| 39 | ARX P2;P1.set_statistics(2, V0, nu0); P2.set_parameters( frg ); |
---|
| 40 | ARX PG;PG.set_statistics(3, V0g, nu0);; PG.set_parameters( frg ); |
---|
[161] | 41 | //Test estimation |
---|
[204] | 42 | int ndat = 200; |
---|
[161] | 43 | int t; |
---|
[162] | 44 | |
---|
[161] | 45 | // Logging |
---|
[213] | 46 | dirfilelog L ( "exp/merg",10 ); |
---|
[161] | 47 | |
---|
[198] | 48 | int Li_Eth1 = L.add ( thri,"P1" ); |
---|
| 49 | int Li_Eth2 = L.add ( thri,"P2" ); |
---|
| 50 | int Li_Ethg = L.add ( thrg,"PG" ); |
---|
| 51 | int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" ); |
---|
| 52 | int Li_LL = L.add ( RV ( "{1 2 G }" ), "LL" ); |
---|
[204] | 53 | int Li_Pred = L.add ( RV ( "{1 2 G ar ge }" ), "Pred" ); |
---|
[162] | 54 | |
---|
| 55 | L.init(); |
---|
| 56 | |
---|
| 57 | vec Yt ( ndat ); |
---|
[213] | 58 | vec yt ( 1 ); |
---|
[162] | 59 | |
---|
| 60 | Yt.set_subvector ( 0,randn ( ord ) ); //initial values |
---|
[270] | 61 | vec rgrg ( thrg._dsize() -1 ); // constant terms are in! |
---|
| 62 | vec rgr1 ( thri._dsize() -1 ); |
---|
| 63 | vec rgr2 ( thri._dsize() -1 ); |
---|
[162] | 64 | |
---|
[213] | 65 | vec PredLLs ( 5 ); |
---|
| 66 | vec PostLLs ( 3 ); |
---|
| 67 | vec PredLLs_m=zeros ( 5 ); |
---|
[198] | 68 | ivec ind_r1 = "0 1 3"; |
---|
| 69 | ivec ind_r2 = "0 2 3"; |
---|
[162] | 70 | for ( t=0; t<ndat; t++ ) { |
---|
[161] | 71 | // True system |
---|
[162] | 72 | if ( t>0 ) { |
---|
| 73 | rgrg ( 0 ) =Yt ( t-1 ); |
---|
| 74 | rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3); |
---|
[213] | 75 | rgrg ( 2 ) = pow(cos ( ( t/60.0 ) *pi ),3); |
---|
| 76 | rgrg ( 3 ) = 1.0; // constant term |
---|
[162] | 77 | |
---|
[213] | 78 | rgr1 ( 0 ) = rgrg ( 0 ); rgr1 ( 1 ) = rgrg ( 1 ); rgr1 ( 2 ) = rgrg ( 3 ); // no u2 |
---|
| 79 | rgr2 ( 0 ) = rgrg ( 0 ); rgr2 ( 1 ) = rgrg ( 2 ); rgr2 ( 2 ) = rgrg ( 3 ); // no u1 |
---|
| 80 | |
---|
[162] | 81 | Yt ( t ) = thg*rgrg + sqr * NorRNG(); |
---|
| 82 | |
---|
[198] | 83 | // Test predictors |
---|
[213] | 84 | if ( t>2 ) { |
---|
[270] | 85 | mlnorm<ldmat>* P1p = P1.predictor();// ( y,concat ( ym,u1 ) ); |
---|
| 86 | mlnorm<ldmat>* P2p = P2.predictor();// ( y,concat ( ym,u2 ) ); |
---|
| 87 | mlnorm<ldmat>* Pgp = PG.predictor();// ( y,concat ( ym,uu ) ); |
---|
[213] | 88 | |
---|
| 89 | Array<mpdf*> A ( 2 ); A ( 0 ) =P1p;A ( 1 ) =P2p; |
---|
| 90 | merger M ( A ); |
---|
[270] | 91 | enorm<ldmat> g0;g0.set_rv ( concat ( yy,uu ) ); g0.set_parameters ( "0 0 0 0 ",3*eye ( 4 ) ); |
---|
[213] | 92 | M.set_parameters ( 1e8, 101,1 ); |
---|
| 93 | M.merge ( &g0 ); |
---|
| 94 | yt ( 0 ) = Yt ( t ); |
---|
| 95 | double P1pl = P1p->evallogcond ( yt,rgr1 ); |
---|
| 96 | double P2pl = P2p->evallogcond ( yt,rgr2 ); |
---|
| 97 | double PGpl = Pgp->evallogcond ( yt,rgrg ); |
---|
[198] | 98 | { |
---|
[213] | 99 | cout << "T: " << t << "yt: " << yt << endl; |
---|
[204] | 100 | cout << "yt_1: " << P1p->_epdf().mean() << endl; |
---|
[213] | 101 | cout << *P1p <<endl; |
---|
[204] | 102 | cout << "yt_2: " << P2p->_epdf().mean() << endl; |
---|
| 103 | cout << "yt_G: " << P2p->_epdf().mean() << endl; |
---|
| 104 | } |
---|
| 105 | double cP1pl; |
---|
| 106 | double cP2pl; |
---|
| 107 | { |
---|
[213] | 108 | ARX* Apred = ( ARX* ) M._Mix()._Coms ( 0 ); |
---|
[270] | 109 | enorm<ldmat>* MP= Apred->epredictor ();//( concat ( yy,uu ) ); |
---|
[213] | 110 | enorm<ldmat>* mP1p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u1 ) ); |
---|
| 111 | enorm<ldmat>* mP2p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u2 ) ); |
---|
| 112 | mlnorm<ldmat>* cP1p = ( mlnorm<ldmat>* ) mP1p->condition ( y ); |
---|
| 113 | mlnorm<ldmat>* cP2p = ( mlnorm<ldmat>* ) mP2p->condition ( y ); |
---|
[198] | 114 | |
---|
[213] | 115 | cP1pl = cP1p->evallogcond ( yt,rgr1 ); |
---|
| 116 | cP2pl = cP2p->evallogcond ( yt,rgr2 ); |
---|
| 117 | |
---|
[204] | 118 | cout << "ytm1: " << cP1p->_epdf().mean() << endl; |
---|
| 119 | cout << "ytm2: " << cP2p->_epdf().mean() << endl; |
---|
[213] | 120 | cout << *cP2p << endl; |
---|
[198] | 121 | } |
---|
| 122 | |
---|
[204] | 123 | PredLLs *=frg; |
---|
[213] | 124 | PredLLs += concat ( vec_3 ( P1pl, P2pl, PGpl ), vec_2 ( cP1pl,cP2pl ) ); |
---|
| 125 | L.logit ( Li_Pred, PredLLs ); //log-normal |
---|
| 126 | |
---|
[198] | 127 | delete P1p; |
---|
| 128 | delete P2p; |
---|
| 129 | delete Pgp; |
---|
| 130 | } |
---|
[162] | 131 | |
---|
[213] | 132 | if ( t<100 ) { |
---|
| 133 | // 1st |
---|
| 134 | P1.bayes ( concat ( Yt ( t ),rgr1 ) ); |
---|
| 135 | // 2nd |
---|
| 136 | P2.bayes ( concat ( Yt ( t ),rgr2 ) ); |
---|
| 137 | |
---|
| 138 | //Global |
---|
| 139 | PG.bayes ( concat ( Yt ( t ),rgrg ) ); |
---|
| 140 | } |
---|
[162] | 141 | //Merger |
---|
| 142 | } |
---|
[271] | 143 | L.logit ( Li_Eth1, P1.posterior().mean() ); |
---|
| 144 | L.logit ( Li_Eth2, P2.posterior().mean() ); |
---|
| 145 | L.logit ( Li_Ethg, PG.posterior().mean() ); |
---|
[198] | 146 | L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) ); |
---|
| 147 | PostLLs *= frg; |
---|
| 148 | PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() ); |
---|
| 149 | L.logit ( Li_LL, PostLLs ); |
---|
[213] | 150 | L.step ( ); |
---|
[161] | 151 | } |
---|
[162] | 152 | L.finalize( ); |
---|
| 153 | L.itsave ( "merg.it" ); |
---|
[161] | 154 | } |
---|