1 | #include <itpp/itmex.h> |
---|
2 | |
---|
3 | #include "../../bdm/estim/arx.h" |
---|
4 | |
---|
5 | using namespace itpp; |
---|
6 | |
---|
7 | void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) { |
---|
8 | // Check the number of inputs and output arguments |
---|
9 | if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" ); |
---|
10 | if ( n_input!=3 ) mexErrMsgTxt ( "Wrong number of input variables!" ); |
---|
11 | |
---|
12 | // Convert input variables to IT++ format |
---|
13 | mat Data = mxArray2mat ( input[0] ); |
---|
14 | if (Data.rows()!=2) |
---|
15 | { if (Data.cols()==2) Data=Data.T(); else mexErrMsgTxt("Data are not 2D!");} |
---|
16 | int ilow = mxArray2int ( input[1] )-1; |
---|
17 | int ihi = mxArray2int ( input[2] )-1; //coreection for different indeces |
---|
18 | |
---|
19 | // ------------------ Start of routine --------------------------- |
---|
20 | // ------------------ Start of routine --------------------------- |
---|
21 | int ndat=Data.cols(); |
---|
22 | int npsi=Data.rows()+1; //add a constant |
---|
23 | RV thr ( "1","{ theta_r }",vec_1 ( npsi ),vec_1 ( 0 ) ); |
---|
24 | |
---|
25 | mat V0=1e-8*eye ( npsi ); V0 ( 0,0 ) =10.0; |
---|
26 | vec rgr(npsi); |
---|
27 | double nu0=2; |
---|
28 | // fitting a linear part => third coef is "1" |
---|
29 | rgr(2) = 1.0; |
---|
30 | |
---|
31 | // RESULTS |
---|
32 | mat Tlls=zeros(ihi,ihi); |
---|
33 | // AR model |
---|
34 | ARX Ar ( thr,V0,nu0 ); |
---|
35 | |
---|
36 | //Initialize |
---|
37 | for (int i=ilow; i<ihi; i++) { |
---|
38 | rgr.set_subvector(0,Data.get_col ( i )); |
---|
39 | Ar.bayes(rgr); |
---|
40 | } |
---|
41 | Ar.get_parameters(V0,nu0); |
---|
42 | V0 /=V0(1,1); |
---|
43 | nu0=3; |
---|
44 | Ar.set_parameters(V0,nu0); |
---|
45 | |
---|
46 | // FLATTEN |
---|
47 | for (int i=ilow; i<ihi; i++) { |
---|
48 | Ar.set_parameters(V0,nu0); |
---|
49 | for (int t=i+1; t<ihi; t++) { |
---|
50 | rgr.set_subvector(0,Data.get_col ( t )); |
---|
51 | Ar.bayes(rgr); |
---|
52 | Tlls(i,t)=Ar._tll(); |
---|
53 | } |
---|
54 | } |
---|
55 | |
---|
56 | max_index(Tlls,ilow,ihi); |
---|
57 | |
---|
58 | // ------------------ End of routine ----------------------------- |
---|
59 | |
---|
60 | output[0] = mxCreateDoubleMatrix ( 2,1, mxREAL ); |
---|
61 | ivec2mxArray (vec_2(ilow+1,ihi+1) , output[0] ); |
---|
62 | |
---|
63 | if ( n_output>1 ) { |
---|
64 | // Create output vectors |
---|
65 | output[1] = mxCreateDoubleMatrix ( Tlls.rows(),Tlls.cols(), mxREAL ); |
---|
66 | // Convert the IT++ format to Matlab format for output |
---|
67 | mat2mxArray ( Tlls, output[1] ); |
---|
68 | } |
---|
69 | |
---|
70 | } |
---|