root/applications/robust/main.cpp @ 1376

Revision 1376, 24.8 kB (checked in by sindj, 13 years ago)

Uprava integrace v souladu s clankem. JS

RevLine 
[976]1
2/*!
3\file
4\brief Robust
5\author Vasek Smidl
6
7 */
8
[1337]9#include "estim/arx.h"
[976]10#include "robustlib.h"
[1216]11#include <vector>
[1284]12#include <iostream>
[1282]13#include <fstream>
[1376]14//#include <itpp/itsignal.h>
[1361]15#include "windows.h"
16#include "ddeml.h"
17#include "stdio.h"
[1282]18
[1361]19//#include "DDEClient.h"
20//#include <conio.h>
[1358]21
[1361]22
[1208]23using namespace itpp;
[1337]24using namespace bdm;
[976]25
[1361]26//const int emlig_size = 2;
27//const int utility_constant = 5;
[1268]28
[1361]29const int max_model_order = 2;
[1376]30const double apriorno=0.02;
[1272]31
[1376]32/*
[1361]33HDDEDATA CALLBACK DdeCallback(
[1376]34        UINT uType,     // Transaction type.
35        UINT uFmt,      // Clipboard data format.
36        HCONV hconv,    // Handle to the conversation.
37        HSZ hsz1,       // Handle to a string.
38        HSZ hsz2,       // Handle to a string.
39        HDDEDATA hdata, // Handle to a global memory object.
40        DWORD dwData1,  // Transaction-specific data.
41        DWORD dwData2)  // Transaction-specific data.
[1361]42{
[1376]43        return 0;
[1361]44}
45
46void DDERequest(DWORD idInst, HCONV hConv, char* szItem)
47{
48        HSZ hszItem = DdeCreateStringHandle(idInst, szItem, 0);
49        HDDEDATA hData = DdeClientTransaction(NULL,0,hConv,hszItem,CF_TEXT,
[1365]50                                                                        XTYP_ADVSTART,TIMEOUT_ASYNC , NULL); //TIMEOUT_ASYNC
[1361]51        if (hData==NULL)
52        {
53                printf("Request failed: %s\n", szItem);
54        }
55
56        if (hData==0)
57        {
58                printf("Request failed: %s\n", szItem);
59        }
60}
61
[1367]62DWORD WINAPI ThrdFunc( LPVOID n )
63{   
64    return 0;
65}
[1376]66*/
[1367]67
[1361]68class model
69{
70public:
[1376]71        set<pair<int,int>> ar_components;
[1358]72
[1361]73        // Best thing would be to inherit the two models from a single souce, this is planned, but now structurally
74        // problematic.
[1376]75        RARX* my_rarx; //vzmenovane parametre pre triedu model
76        ARX* my_arx;
[1361]77
78        bool has_constant;
[1376]79        int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit
[1361]80        int  predicted_channel;
81        mat* data_matrix;
82       
[1376]83        model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor
[1361]84                  bool robust, 
85                  bool has_constant, 
86                  int window_size, 
[1376]87                  int predicted_channel,
[1361]88                  mat* data_matrix)
[1358]89        {
[1376]90                this->ar_components.insert(ar_components.begin(),ar_components.end());
91                this->has_constant      = has_constant;
92                this->window_size       = window_size;
93                this->predicted_channel = predicted_channel;
94                this->data_matrix       = data_matrix;
[1361]95
96                if(robust)
97                {
98                        if(has_constant)
99                        {
100                                my_rarx = new RARX(ar_components.size()+1,window_size,true);
101                                my_arx  = NULL;
102                        }
[1376]103                else
[1361]104                        {
105                                my_rarx = new RARX(ar_components.size(),window_size,false);
106                                my_arx  = NULL;
107                        }
108                }
109                else
110                {
111                        my_rarx = NULL;
[1376]112                        my_arx  = new ARX();
[1361]113                        mat V0;
114
115                        if(has_constant)
116                        {                               
[1376]117                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst
[1370]118                                V0(0,0) = 1;
[1361]119                                my_arx->set_constant(true);     
120                               
121                        }
122                        else
123                        {
124                               
[1376]125                                V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
[1370]126                                V0(0,0) = 1;
[1361]127                                my_arx->set_constant(false);                           
128                               
129                        }
130
131                        my_arx->set_statistics(1, V0); 
[1376]132                        //my_arx->set_parameters(window_size);
[1361]133                        my_arx->validate();
134                }
[1358]135        }
[1361]136
[1376]137        void data_update(int time) //vlozime cas a ono vlozi do data_vector podmineky(conditions) a predikce, ktore pouzije do bayes
[1358]138        {
[1376]139                vec data_vector;
140                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
141                {  //ar?iterator ide len od 1 pod 2, alebo niekedy len 1
142                        data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
143// do data vector vlozi pre dany typ regresoru prislusne cisla z data_matrix. Ale ako? preco time-ar_iterator->second?
[1358]144                }
[1376]145                if(my_rarx!=NULL)
146                {       //pre robusr priradi az tu do data_vector aj rpedikciu
147                        data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
148                        my_rarx->bayes(data_vector);
149                }
[1358]150                else
151                {
[1376]152                        vec pred_vec;//tu sa predikcia zadava zvlast
153                        pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
154                        my_arx->bayes(pred_vec,data_vector);
[1361]155                }
156        }
157
[1376]158        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to
[1367]159        {
[1376]160                vec condition_vector;
161                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
[1367]162                {
[1376]163                        condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1));
164                }
[1367]165
[1376]166                if(my_rarx!=NULL)
167                {
168                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(sample_size);
[1367]169
[1376]170                        //cout << imp_samples.first << endl;
171                       
172                        vec sample_prediction;                 
173                        for(int t = 0;t<sample_size;t++)
[1367]174                        {
[1376]175                                vec lap_sample = condition_vector;
[1367]176                               
[1376]177                                if(has_constant)
[1367]178                                {
[1376]179                                        lap_sample.ins(lap_sample.size(),1.0);
[1367]180                                }
[1376]181                               
182                                lap_sample.ins(0,(*LapRNG)());
[1367]183
[1376]184                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                             
[1367]185                        }
186
[1376]187                        return pair<vec,vec>(imp_samples.first,sample_prediction);
188                }
189                else
190                {
191                        mat samples = my_arx->posterior().sample_mat(sample_size);
192                       
193                        vec sample_prediction;
194                        for(int t = 0;t<sample_size;t++)
195                        {
196                                vec gau_sample = condition_vector;
[1367]197                               
[1376]198                                if(has_constant)
[1367]199                                {
[1376]200                                        gau_sample.ins(gau_sample.size(),1.0);
[1367]201                                }
[1376]202                               
203                                gau_sample.ins(0,randn());
[1367]204
[1376]205                                sample_prediction.ins(0,gau_sample*samples.get_col(t));                         
[1367]206                        }
[1376]207
208                        return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction);
[1367]209                }
210       
211        }
212
213
[1376]214        static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels)
[1361]215        {
[1376]216                set<set<pair<int,int>>> created_model_types;           
[1361]217
[1376]218                if(max_order == 1)//ukoncovacia vetva
[1361]219                {                       
[1376]220                        for(int channel = 0;channel<number_of_channels;channel++)//pre AR 1 model vytvori kombinace kanalov v prvom kroku poyadu
[1358]221                        {
[1376]222                                set<pair<int,int>> returned_type;
223                                returned_type.insert(pair<int,int>(channel,1));  //??
224                                created_model_types.insert(returned_type);
[1358]225                        }
[1361]226
227                        return created_model_types;
228                }
229                else
230                {
[1376]231                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);//tu uz mame ulozene kombinace o jeden krok dozadu  //rekuryivne volanie
232                        set<set<pair<int,int>>> returned_types;
[1361]233                       
[1376]234                        for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++)
[1361]235                        {                               
236                               
237                                for(int order = 1; order<=max_order; order++)
[1358]238                                {
[1361]239                                        for(int channel = 0;channel<number_of_channels;channel++)
240                                        {
[1376]241                                                set<pair<int,int>> returned_type;
242                                                pair<int,int> new_pair = pair<int,int>(channel,order);//??
243                                                if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) //??
[1361]244                                                {
[1376]245                                                        returned_type.insert((*model_ref).begin(),(*model_ref).end()); //co vlozi na zaciatok retuned_type?
246                                                        returned_type.insert(new_pair);
247                                                       
248
249                                                        returned_types.insert(returned_type);                                                   
[1361]250                                                }
251                                        }
[1358]252                                }
253                        }
[1361]254
[1376]255                        created_model_types.insert(returned_types.begin(),returned_types.end());
[1361]256
257                        return created_model_types;
258                }               
[1358]259        }
[1361]260};
261
262
263
264
[1367]265int main ( int argc, char* argv[] ) {   
[1358]266       
[1367]267        /*
[1376]268        DWORD Id;
[1367]269        HANDLE hThrd = CreateThread( NULL, 0, (LPTHREAD_START_ROUTINE)ThrdFunc, (LPVOID)1, 0, &Id);
[1376]270 
271        if ( !hThrd )
[1367]272    {
273        cout<<"Error Creating Threads,,,,.exiting"<<endl;
274        return -1;
275    }
276    Sleep ( 100 );
[1376]277
[1361]278       
[1376]279        char szApp[] = "MT4";
280        char szTopic[] = "QUOTE";       
281        char szItem1[] = "EURUSD";     
[1361]282
283        //DDE Initialization
[1376]284        DWORD idInst=0;
[1361]285        UINT iReturn;
286        iReturn = DdeInitialize(&idInst, (PFNCALLBACK)DdeCallback,
287                                                        APPCLASS_STANDARD | APPCMD_CLIENTONLY, 0 );
288        if (iReturn!=DMLERR_NO_ERROR)
289        {
290                printf("DDE Initialization Failed: 0x%04x\n", iReturn);
291                Sleep(1500);
292                return 0;
293        }
294
295        //DDE Connect to Server using given AppName and topic.
296        HSZ hszApp, hszTopic;
297        HCONV hConv;
298        hszApp = DdeCreateStringHandle(idInst, szApp, 0);
299        hszTopic = DdeCreateStringHandle(idInst, szTopic, 0);
300        hConv = DdeConnect(idInst, hszApp, hszTopic, NULL);
[1367]301        //DdeFreeStringHandle(idInst, hszApp);
302        //DdeFreeStringHandle(idInst, hszTopic);
[1361]303        if (hConv == NULL)
304        {
305                printf("DDE Connection Failed.\n");
306                Sleep(1500); DdeUninitialize(idInst);
307                return 0;
308        }
309
310        //Execute commands/requests specific to the DDE Server.
311       
[1376]312        DDERequest(idInst, hConv, szItem1);             
[1361]313       
[1367]314        while(1)
315        {
[1376]316                MSG    msg;
317                BOOL   MsgReturn = GetMessage ( &msg , NULL , 0 , 0 );
318           
319                if(MsgReturn)
[1367]320                {
[1376]321                        TranslateMessage(&msg);
322                        DispatchMessage(&msg);                 
[1367]323                }
324        }
[1376]325
[1361]326        //DDE Disconnect and Uninitialize.
[1367]327        DdeDisconnect(hConv);
328        DdeUninitialize(idInst);
[1376]329        */
[1361]330
[1376]331       
[1361]332
[1376]333        /*
334        // EXPERIMENT: 100 AR model generated time series of length of 30 from y_t=0.95*y_(t-1)+0.05*y_(t-2)+0.2*e_t,
335        // where e_t is normally, student(4) and cauchy distributed are tested using robust AR model, to obtain the
336        // variance of location parameter estimators and compare it to the classical setup.
337        vector<vector<vector<string>>> string_lists;
338        string_lists.push_back(vector<vector<string>>());
339        string_lists.push_back(vector<vector<string>>());
340        string_lists.push_back(vector<vector<string>>());
[1186]341
[1376]342        char* file_strings[3] = {"c:\\ar_normal.txt", "c:\\ar_student.txt", "c:\\ar_cauchy.txt"};
343       
[1282]344
[1376]345        for(int i = 0;i<3;i++)
346        {       
347                ifstream myfile(file_strings[i]);
[1282]348                if (myfile.is_open())
[1376]349                {
350                        while ( myfile.good() )
351                        {
352                                string line;
353                                getline(myfile,line);
354                               
355                                vector<string> parsed_line;
[1282]356                                while(line.find(',') != string::npos)
357                                {
[1376]358                                        int loc = line.find(',');
359                                        parsed_line.push_back(line.substr(0,loc));
[1282]360                                        line.erase(0,loc+1);                                   
[1376]361                                }                               
[1282]362
[1376]363                                string_lists[i].push_back(parsed_line);
364                        }
365                        myfile.close();
[1282]366                }
[1376]367        }
[1282]368
[1376]369        for(int j = 0;j<string_lists.size();j++)
370        {
[1301]371               
[1376]372                for(int i = 0;i<string_lists[j].size()-1;i++)
373                {
374                        vector<vec> conditions;
375                        //emlig* emliga = new emlig(2);
376                        RARX* my_rarx = new RARX(2,30);
[1301]377
[1376]378                        for(int k = 1;k<string_lists[j][i].size();k++)
[1282]379                        {
[1376]380                                vec condition;
381                                //condition.ins(0,1);                           
382                                condition.ins(0,string_lists[j][i][k]);                         
383                                conditions.push_back(condition);
[1282]384
[1376]385                                //cout << "orig:" << condition << endl;
386
387                                if(conditions.size()>1)
388                                {               
389                                        conditions[k-2].ins(0,string_lists[j][i][k]);
390                                       
[1282]391                                }
[1376]392
393                                if(conditions.size()>2)
394                                {
395                                        conditions[k-3].ins(0,string_lists[j][i][k]);
396
397                                        //cout << "modi:" << conditions[k-3] << endl;
398
399                                        my_rarx->bayes(conditions[k-3]);
400
401                                       
402                                        //if(k>5)
403                                        //{
404                                        //      cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
405                                        //}
406                                       
407                                }                               
408                               
[1282]409                        }
410
[1376]411                        //emliga->step_me(0);
412                        /*
413                        ofstream myfile;
414                        myfile.open("c:\\robust_ar1.txt",ios::app);
415                        myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";";
416                        myfile.close();
417
418                        myfile.open("c:\\robust_ar2.txt",ios::app);
419                        myfile << emliga->minimal_vertex->get_coordinates()[1] << ";";
420                        myfile.close();
[1301]421                       
[1284]422
[1376]423                        cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
424                        cout << "Step: " << i << endl;
425                }
[1282]426
[1376]427                cout << "One experiment finished." << endl;
[1284]428
[1376]429                ofstream myfile;
430                myfile.open("c:\\robust_ar1.txt",ios::app);
431                myfile << endl;
432                myfile.close();
[1284]433
[1376]434                myfile.open("c:\\robust_ar2.txt",ios::app);
435                myfile << endl;
436                myfile.close();
437        }*/   
438       
439        // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from
440    // y_t=0.95*y_(t-1)+0.05*y_(t-2)+0.2*e_t, where e_t is normally, student(4) and cauchy distributed. It
441    // can be compared to the classical setup.
442       
443        itpp::Laplace_RNG LapRNG = Laplace_RNG();       
[1301]444
[1376]445        vector<vector<string>> strings;
[1301]446
[1376]447        char* file_string = "c:\\ar_cauchy_single"; 
[1301]448
[1376]449        char dfstring[80];
450        strcpy(dfstring,file_string);
451        strcat(dfstring,".txt");
452       
453       
454        mat data_matrix;
455        ifstream myfile(dfstring);
456        if (myfile.is_open())
457        {               
458                string line;
459                while(getline(myfile,line))
460                {                       
461                        vec data_vector;
462                        while(line.find(',') != string::npos) //zmenil som ciarku za medzeru
463                        {
464                                line.erase(0,1); // toto som sem pridal
465                                int loc2 = line.find('\n');
466                                int loc  = line.find(',');
467                                data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                           
468                                line.erase(0,loc+1);                                   
469                        }
[1301]470
[1376]471                        data_matrix.ins_row(data_matrix.rows(),data_vector);
472                }               
[1361]473
[1376]474                myfile.close(); 
475        }
476        else
477        {
478                cout << "Can't open data file!" << endl;
479        }
480       
481        //konec nacitavania dat
482        set<set<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,data_matrix.rows()); //volanie funkce kde robi kombinace modelov
483                                                                                                //to priradime do model_types, data_matrix.row urcuje pocet kanalov dat
484        vector<model*> models;
485        for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++)
486        {// prechadza rozne typy kanalov, a poctu regresorov
487                for(int window_size = 100;window_size < 101;window_size++)
488                {
489                        models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy
490                        models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));
491                        models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix));
492                        models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));             
493                }
[1365]494
[1376]495                set<pair<int,int>> empty_list;
496                models.push_back(new model(empty_list,false,true,100,0,&data_matrix));
497        }
498       
499        mat result_lognc;
500        // mat result_preds;
501
502        for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols()
503        {       //pocet stlpcov data_matrix je pocet casovych krokov
504                vec cur_res_lognc;
505                // vec preds;
506                vector<string> nazvy;
507                for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++)
508                {//posuvam s apo models, co je pole modelov urobene o cyklus vyssie. Teda som v case time a robim to tam pre vsetky typy modelov, kombinace regresorov
509                        (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne
510                        //if (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu?
511
512                        if((*model_ref)->my_rarx!=NULL) //vklada normalizacnz faktor do cur_res_lognc
513                        {
514                                cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_rarx->posterior->log_nc);
[1361]515                        }
[1376]516                        else
517                        {
518                                cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_arx->posterior().lognc());
519                        }
[1367]520
[1376]521                        // pair<vec,vec> predictions = (*model_ref)->predict(200,time,&LapRNG);
[1367]522
[1376]523                        // preds.ins(preds.size(),(predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size())));
524                        // preds.ins(0,data_matrix.get(0,time+1));
[1361]525
[1376]526                }
[1367]527
[1376]528                result_lognc.ins_col(result_lognc.cols(),cur_res_lognc);
529                // result_preds.ins_col(result_preds.cols(),preds);
530
531                // cout << "Updated." << endl;
532       
[1361]533                /*
[1301]534                vector<vec> conditions;
535                //emlig* emliga = new emlig(2);
[1358]536                RARX* my_rarx = new RARX(2,10,false);
[1337]537               
[1338]538               
[1337]539                mat V0 = 0.0001 * eye ( 3 );
[1349]540                ARX* my_arx = new ARX(0.85);
[1337]541                my_arx->set_statistics ( 1, V0 ); //nu is default (set to have finite moments)
542                my_arx->set_constant ( false );
543                my_arx->validate();
[1338]544               
[1301]545
[1338]546                for(int k = 1;k<strings[j].size();k++)
[1301]547                {
548                        vec condition;
549                        //condition.ins(0,1);                           
550                        condition.ins(0,strings[j][k]);                         
551                        conditions.push_back(condition);
552
553                        //cout << "orig:" << condition << endl;
554
555                        if(conditions.size()>1)
556                        {               
557                                conditions[k-2].ins(0,strings[j][k]);
558                                       
559                        }
560
561                        if(conditions.size()>2)
562                        {
563                                conditions[k-3].ins(0,strings[j][k]);
564
[1349]565                                // cout << "Condition:" << conditions[k-3] << endl;
[1301]566
567                                my_rarx->bayes(conditions[k-3]);
[1338]568                                //my_rarx->posterior->step_me(1);
[1337]569                               
570                                vec cond_vec;
571                                cond_vec.ins(0,conditions[k-3][0]);
572                               
[1338]573                                my_arx->bayes(cond_vec,conditions[k-3].right(2));
[1301]574                                       
[1361]575                                /*
[1346]576                                if(k>8)
[1301]577                                {
[1324]578                                        //my_rarx->posterior->step_me(0);
579
[1346]580                                        //mat samples = my_rarx->posterior->sample_mat(10);
[1343]581
[1346]582                                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(1000);
[1343]583
[1346]584                                        //cout << imp_samples.first << endl;
[1336]585                                       
[1337]586                                        vec sample_prediction;
[1358]587                                        vec averaged_params = zeros(imp_samples.second.rows());
[1346]588                                        for(int t = 0;t<imp_samples.first.size();t++)
[1337]589                                        {
590                                                vec lap_sample = conditions[k-3].left(2);
[1346]591                                                //lap_sample.ins(lap_sample.size(),1.0);
[1337]592                                               
593                                                lap_sample.ins(0,LapRNG());
594
[1346]595                                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));
[1358]596
597                                                averaged_params += imp_samples.first[t]*imp_samples.second.get_col(t);
[1337]598                                        }
599
[1358]600                                        averaged_params = averaged_params*(1/(imp_samples.first*ones(imp_samples.first.size())));
601
602                                        // cout << "Averaged estimated parameters: " << averaged_params << endl;
[1338]603                                       
[1358]604                                        vec sample_pow = sample_prediction;                                     
[1343]605                                       
606                                        // cout << sample_prediction << endl;
[1337]607                                        vec poly_coefs;
[1346]608                                        double prediction;
[1337]609                                        bool stop_iteration = false;
[1343]610                                        int en = 1;
[1337]611                                        do
612                                        {
[1346]613                                                double poly_coef = imp_samples.first*sample_pow/(imp_samples.first*ones(imp_samples.first.size()));
[1337]614
[1346]615                                                if(en==1)
616                                                {
617                                                        prediction = poly_coef;
618                                                }
619
[1343]620                                                poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2);
621
[1337]622                                                if(abs(poly_coef)>numeric_limits<double>::epsilon())
623                                                {
624                                                        sample_pow = elem_mult(sample_pow,sample_prediction);
[1343]625                                                        poly_coefs.ins(0,pow(-1.0,en+1)*poly_coef);
[1337]626                                                }
627                                                else
628                                                {
629                                                        stop_iteration = true;
630                                                }
631                                               
632                                                en++;
[1343]633
634                                                if(en>20)
635                                                {
636                                                        stop_iteration = true;
637                                                }
[1337]638                                        }
639                                        while(!stop_iteration);
640
[1343]641                                        /*
642                                        ofstream myfile_coef;                                           
643
644                                        myfile_coef.open("c:\\coefs.txt",ios::app);
645                                       
646                                        for(int t = 0;t<poly_coefs.size();t++)
647                                        {
648                                                myfile_coef << poly_coefs[t] << ",";                                   
649                                        }
650
651                                        myfile_coef << endl;
652                                        myfile_coef.close();
653                                        */
654
[1349]655                                        //cout << "Coefficients: " << poly_coefs << endl;
[1338]656                                                                               
[1343]657                                        /*
658                                        vec bas_coef = vec("1.0 2.0 -8.0");
659                                        cout << "Coefs: " << bas_coef << endl;
660                                        cvec actions2 = roots(bas_coef);
661                                        cout << "Roots: " << actions2 << endl;
662                                        */
663                                       
[1361]664                                    /*
[1346]665
[1338]666                                        cvec actions = roots(poly_coefs);
[1343]667                                       
668
[1338]669                                        bool is_max = false;
670                                        for(int t = 0;t<actions.size();t++)
671                                        {
[1343]672                                                if(actions[t].imag() == 0)
[1338]673                                                {
[1343]674                                                        double second_derivative = 0;
675                                                        for(int q = 1;q<poly_coefs.size();q++)
676                                                        {
677                                                                second_derivative+=poly_coefs[q]*pow(actions[t].real(),q-1)*q;
678                                                        }
679
680                                                        if(second_derivative<0)
681                                                        {
682                                                                cout << "Action:" << actions[t].real() << endl;
683
684                                                                is_max = true;
685                                                        }
[1338]686                                                }
687                                        }
[1301]688
[1338]689                                        if(!is_max)
690                                        {
691                                                cout << "No maximum." << endl;
692                                        }
693
694                                        // cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
695
[1346]696                                        /*
[1337]697                                        double prediction = 0;
698                                        for(int s = 1;s<samples.rows();s++)
[1336]699                                        {
700                                               
[1346]701                                                double avg_parameter = imp_samples.get_row(s)*ones(samples.cols())/samples.cols();
[1337]702
703                                                prediction += avg_parameter*conditions[k-3][s-1];
704
[1336]705                                               
[1337]706                                               
707                                                /*
[1336]708                                                ofstream myfile;
709                                                char fstring[80];
710                                                strcpy(fstring,file_strings[j]);
[1301]711
[1336]712                                                char es[5];
713                                                strcat(fstring,itoa(s,es,10));
714
715                                                strcat(fstring,"_res.txt");
716                                               
717
718                                                myfile.open(fstring,ios::app);
719                                               
720                                                //myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
721                                                myfile << avg_parameter;
722                                               
723                                                if(k!=strings[j].size()-1)
724                                                {
725                                                        myfile << ",";
726                                                }
727                                                else
728                                                {
729                                                        myfile << endl;
730                                                }
731                                                myfile.close();
[1337]732                                                */
733
[1338]734                                       
[1346]735                                        //}
736
737                                        // cout << "Prediction: "<< prediction << endl;
[1361]738                                        /*
[1337]739                                        enorm<ldmat>* pred_mat = my_arx->epredictor(conditions[k-3].left(2));
740                                        double prediction2 = pred_mat->mean()[0];
[1361]741                                        */
[1337]742
[1365]743                                       
[1376]744                                        ofstream myfile;
745                                        char fstring[80];                                       
746                                        strcpy(fstring,file_string);
[1337]747                                       
[1376]748                                        strcat(fstring,"lognc.txt");                                   
749
750                                        myfile.open(fstring,ios::app);
751                                       
752                                        // myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
753                                       
754                                        if(time == max_model_order)
755                                        { 
756                                                for(int i = 0;i<cur_res_lognc.size();i++)
757                                                {
758                                                        for(set<pair<int,int>>::iterator ar_ref = models[i]->ar_components.begin();ar_ref != models[i]->ar_components.end();ar_ref++)
759                                                        {
760                                                                myfile << (*ar_ref).second << (*ar_ref).first;                                                 
761                                                        }
762
763                                                        myfile << ".";
764
765                                                        if(models[i]->my_arx == NULL)
766                                                        {
767                                                                myfile << "1";
768                                                        }
769                                                        else
770                                                        {
771                                                                myfile << "0";
772                                                        }
773                                                               
774                                                        if(models[i]->has_constant)
775                                                        {
776                                                                myfile << "1";
777                                                        }
778                                                        else
779                                                        {
780                                                                myfile << "0";
781                                                        }
782
783                                                        myfile << ",";
784                                                }
785
786                                                myfile << endl;
787                                        }
788                                       
789                                        for(int i = 0;i<cur_res_lognc.size();i++)
790                                        {
791                                                myfile << cur_res_lognc[i] << ' ';//zmenil som ciarku ze medzeru
792                                        }
793
794                                        myfile << endl;                         
795                                       
796                                        myfile.close();
797                        }
[1361]798                                        /*
[1337]799                                        myfile.open(f2string,ios::app);
800                                        myfile << prediction2;
801                                       
802                                        if(k!=strings[j].size()-1)
803                                        {
804                                                myfile << ",";
805                                        }
806                                        else
807                                        {
808                                                myfile << endl;
809                                        }
810                                        myfile.close();
[1361]811                                        //*//*
[1337]812
[1319]813                                }                                       
[1361]814                        }       */
[1301]815                       
816                        //emliga->step_me(0);
817                        /*
818                        ofstream myfile;
819                        myfile.open("c:\\robust_ar1.txt",ios::app);
820                        myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";";
821                        myfile.close();
822
823                        myfile.open("c:\\robust_ar2.txt",ios::app);
824                        myfile << emliga->minimal_vertex->get_coordinates()[1] << ";";
825                        myfile.close();
826                       
827
828                        cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
829                        cout << "Step: " << i << endl;*/
[1361]830                //}
[1337]831
832
[1365]833        //}
[1337]834
835
836        // EXPERIMENT: One step ahead price prediction. Comparison of classical and robust model using optimal trading
837    //             with maximization of logarithm of one-step ahead wealth.
838
839       
[1301]840               
841                /*
842                cout << "One experiment finished." << endl;
843
844                ofstream myfile;
845                myfile.open("c:\\robust_ar1.txt",ios::app);
846                myfile << endl;
847                myfile.close();
848
849                myfile.open("c:\\robust_ar2.txt",ios::app);
850                myfile << endl;
851                myfile.close();*/
[1300]852       
[1301]853
854        //emlig* emlig1 = new emlig(emlig_size);
855
856        //emlig1->step_me(0);
857        //emlig* emlig2 = new emlig(emlig_size);
[1300]858       
[1267]859        /*
860        emlig1->set_correction_factors(4);
[1266]861
[1267]862        for(int j = 0;j<emlig1->correction_factors.size();j++)
863        {
864                for(set<my_ivec>::iterator vec_ref = emlig1->correction_factors[j].begin();vec_ref!=emlig1->correction_factors[j].end();vec_ref++)
865                {
[1268]866                        cout << j << "    ";
867                       
[1267]868                        for(int i=0;i<(*vec_ref).size();i++)
869                        {
870                                cout << (*vec_ref)[i];
871                        }
872
873                        cout << endl;
874                }
[1268]875        }*/
876       
[1301]877        /*
[1300]878    vec condition5 = "1.0 1.0 1.01";//"-0.3 1.7 1.5";
879
[1299]880        emlig1->add_condition(condition5);
[1301]881        //emlig1->step_me(0);
882
883
884        vec condition1a = "-1.0 1.02 0.5";
[1300]885        //vec condition1b = "1.0 1.0 1.01";
[1301]886        emlig1->add_condition(condition1a);
[1300]887        //emlig2->add_condition(condition1b);
[1267]888
[1301]889        vec condition2a = "-0.3 1.7 1.5";
[1300]890        //vec condition2b = "-1.0 1.0 1.0";
[1301]891        emlig1->add_condition(condition2a);
[1300]892        //emlig2->add_condition(condition2b);
[1234]893
[1301]894        vec condition3a = "0.5 -1.01 1.0";
[1300]895        //vec condition3b = "0.5 -1.01 1.0";
[1280]896
[1301]897        emlig1->add_condition(condition3a);
[1300]898        //emlig2->add_condition(condition3b);   
[1280]899
[1301]900        vec condition4a = "-0.5 -1.0 1.0";
[1300]901        //vec condition4b = "-0.5 -1.0 1.0";   
902
[1301]903        emlig1->add_condition(condition4a);
[1300]904        //cout << "************************************************" << endl;
905        //emlig2->add_condition(condition4b);
906        //cout << "************************************************" << endl;
907       
[1299]908        //cout << emlig1->minimal_vertex->get_coordinates();
[1300]909       
[1301]910        //emlig1->remove_condition(condition3a);
911        //emlig1->step_me(0);
912        //emlig1->remove_condition(condition2a);
913        //emlig1->remove_condition(condition1a);
914        //emlig1->remove_condition(condition5);
915       
[1275]916
[1299]917        //emlig1->step_me(0);
918        //emlig2->step_me(0);
919       
[1282]920
921        // DA SE POUZIT PRO VYPIS DO SOUBORU
[1275]922        // emlig1->step_me(0);
923
924        //emlig1->remove_condition(condition1);
925       
[1301]926       
[1275]927
928       
[1301]929       
[1275]930        /*
[1282]931        for(int i = 0;i<100;i++)
[1219]932        {
[1282]933                cout << endl << "Step:" << i << endl;           
[1208]934
[1268]935                double condition[emlig_size+1];         
[1220]936
[1268]937                for(int k = 0;k<=emlig_size;k++)
938                {
[1272]939                        condition[k] = (rand()-RAND_MAX/2)/1000.0;             
[1268]940                }
941                       
[1216]942
[1268]943                vec* condition_vec = new vec(condition,emlig_size+1);
[1219]944                emlig1->add_condition(*condition_vec);
[1271]945
[1272]946                /*
947                for(polyhedron* toprow_ref = emlig1->statistic.rows[emlig_size]; toprow_ref != emlig1->statistic.end_poly; toprow_ref = toprow_ref->next_poly)
948                {
949                        cout << ((toprow*)toprow_ref)->probability << endl;
950                }
951                */
[1275]952                /*
[1271]953                cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
[1268]954       
[1272]955                /*
[1271]956                if(i-emlig1->number_of_parameters >= 0)
957                {
958                        pause(30);
959                }
[1272]960                */
[1219]961
[1271]962                // emlig1->step_me(i);
[1219]963               
[1272]964                /*
[1219]965                vector<int> sizevector;
966                for(int s = 0;s<=emlig1->number_of_parameters;s++)
967                {
968                        sizevector.push_back(emlig1->statistic_rowsize(s));
969                }
[1272]970                */
[1275]971        //}
972   
[1219]973
974
975       
976        /*
977        emlig1->step_me(1);
978
979        vec condition = "2.0 0.0 1.0"; 
980
[1208]981        emlig1->add_condition(condition);
982
[1216]983        vector<int> sizevector;
984        for(int s = 0;s<=emlig1->number_of_parameters;s++)
985        {
986                sizevector.push_back(emlig1->statistic_rowsize(s));
987        }
988
[1219]989        emlig1->step_me(2);
[1216]990
[1219]991        condition = "2.0 1.0 0.0";
[1216]992
993        emlig1->add_condition(condition);
994
995        sizevector.clear();
996        for(int s = 0;s<=emlig1->number_of_parameters;s++)
997        {
998                sizevector.push_back(emlig1->statistic_rowsize(s));
999        }
[1219]1000        */
[1216]1001
[1376]1002        return 0;
1003}
[976]1004
[1282]1005
Note: See TracBrowser for help on using the browser.