root/applications/robust/main.cpp @ 1412

Revision 1412, 18.3 kB (checked in by sindj, 12 years ago)

Main upravena pro experimenty na burzovnich datech. JS

Line 
1
2/*!
3\file
4\brief Robust
5\author Vasek Smidl
6
7 */
8
9#include "estim/arx.h"
10#include "robustlib.h"
11#include <vector>
12#include <iostream>
13#include <fstream>
14//#include <itpp/itsignal.h>
15#include "windows.h"
16#include "ddeml.h"
17#include "stdio.h"
18#include <dirent.h>
19//#include <itpp/itoptim.h>
20
21//#include "DDEClient.h"
22//#include <conio.h>
23
24
25using namespace itpp;
26using namespace bdm;
27
28//const int emlig_size = 2;
29//const int utility_constant = 5;
30
31const int max_model_order         = 2;
32const double apriorno             = 0.01;
33const int max_window_size         = 40;
34const int utility_order           = 25;
35const int prediction_time         = 30;
36const double min_utility_argument = 0.001;
37const double max_investment       = 10.0;
38const char* commodity                 = "CL\\";
39
40/*
41HDDEDATA CALLBACK DdeCallback(
42        UINT uType,     // Transaction type.
43        UINT uFmt,      // Clipboard data format.
44        HCONV hconv,    // Handle to the conversation.
45        HSZ hsz1,       // Handle to a string.
46        HSZ hsz2,       // Handle to a string.
47        HDDEDATA hdata, // Handle to a global memory object.
48        DWORD dwData1,  // Transaction-specific data.
49        DWORD dwData2)  // Transaction-specific data.
50{
51        return 0;
52}
53
54void DDERequest(DWORD idInst, HCONV hConv, char* szItem)
55{
56        HSZ hszItem = DdeCreateStringHandle(idInst, szItem, 0);
57        HDDEDATA hData = DdeClientTransaction(NULL,0,hConv,hszItem,CF_TEXT,
58                                                                        XTYP_ADVSTART,TIMEOUT_ASYNC , NULL); //TIMEOUT_ASYNC
59        if (hData==NULL)
60        {
61                printf("Request failed: %s\n", szItem);
62        }
63
64        if (hData==0)
65        {
66                printf("Request failed: %s\n", szItem);
67        }
68}
69
70DWORD WINAPI ThrdFunc( LPVOID n )
71{   
72    return 0;
73}
74*/
75vector<char*> listFiles(char* dir){
76        vector<char*> files;
77                DIR *pDIR;
78        struct dirent *entry;
79        if( pDIR=opendir(dir) ){
80                while(entry = readdir(pDIR)){
81                        if( strcmp(entry->d_name, ".") != 0 && strcmp(entry->d_name, "..") != 0 )
82                                                {                                                       
83                                                        char *file_string = new char[strlen(entry->d_name) + 1];
84                                                        strcpy(file_string, entry->d_name);                                                                                             
85                                                        files.push_back(file_string);
86                                                }
87                }
88                closedir(pDIR);
89        }
90                return files;
91}
92
93double valueCRRAUtility(const double &position, const vec &samples, const int order)
94{
95        double value = 0;
96                       
97        for(int i=0;i<samples.length();i++)
98        {
99                double sample = samples.get(i);
100                if((position*sample+1)>min_utility_argument)
101                {
102                        value += sample/pow(position*sample+1,order+1);
103                }
104                else
105                {
106                        value += (min_utility_argument-1)/position/pow(min_utility_argument,order+1);
107                }
108        }       
109       
110        return value;
111}
112
113double gradientCRRAUtility(const double &position, const vec &samples, const int order)
114{       
115        double value = 0;
116       
117        for(int i=0;i<samples.length();i++)
118        {
119                double sample = samples.get(i);
120                if((position*sample+1)>min_utility_argument)
121                {
122                        value += (-(order+1)*pow(sample,2))/pow(position*sample+1,order+2);
123                }
124        }
125
126        return value;
127}
128
129double newtonRaphson(double startingPoint, double epsilon, vec samples, int order)
130{
131        if(samples.length()>800)
132        {
133                samples.del(801,samples.size()-1);
134        }
135       
136        int count = 0;
137       
138        bool epsilon_reached = false;
139
140        while(count<1000&&!epsilon_reached)
141        {
142                double cur_value = valueCRRAUtility(startingPoint,samples,order);
143                double cur_gradient = gradientCRRAUtility(startingPoint,samples,order);
144
145                startingPoint = startingPoint - cur_value/cur_gradient;
146
147                if(cur_value<epsilon)
148                {
149                        epsilon_reached = true;
150                }
151        }
152
153        if(count==100)
154        {
155                return startingPoint; // can be different!
156        }
157        else
158        {
159                return startingPoint;
160        }
161}
162
163class model
164{
165public:
166        set<pair<int,int>> ar_components;
167
168        // Best thing would be to inherit the two models from a single souce, this is planned, but now structurally
169        // problematic.
170        RARX* my_rarx; //vzmenovane parametre pre triedu model
171        ARXwin* my_arx;
172
173        bool has_constant;
174        int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit
175        int  predicted_channel;
176        mat* data_matrix;
177        vec  predictions;
178        char name[80];
179
180        double previous_lognc;
181       
182        model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor
183                  bool robust, 
184                  bool has_constant, 
185                  int window_size, 
186                  int predicted_channel,
187                  mat* data_matrix)
188        {               
189                this->ar_components.insert(ar_components.begin(),ar_components.end());
190
191                strcpy(name,"M");
192
193                for(set<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++)
194                {
195                        char buffer1[2];
196                        char buffer2[2];
197                        itoa((*ar_ref).first,buffer1,10);
198                        itoa((*ar_ref).second,buffer2,10);
199
200                        strcat(name,buffer1);
201                        strcat(name,buffer2);
202                        strcat(name,"_");
203                }
204
205                this->has_constant      = has_constant;
206
207                if(has_constant)
208                {
209                        strcat(name,"C");
210                }
211
212                this->window_size       = window_size;
213                this->predicted_channel = predicted_channel;
214                this->data_matrix       = data_matrix;
215
216                if(robust)
217                {
218                        previous_lognc = 0;
219                        strcat(name,"R");
220
221                        if(has_constant)
222                        {
223                                my_rarx = new RARX(ar_components.size()+1,window_size,true,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+3);
224                                my_arx  = NULL;
225                        }
226                else
227                        {
228                                my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+2);
229                                my_arx  = NULL;
230                        }
231                }
232                else
233                {
234                        my_rarx = NULL;
235                        my_arx  = new ARXwin();
236                        mat V0;
237
238                        if(has_constant)
239                        {                               
240                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst
241                                V0(0,0) = 0;
242                                my_arx->set_constant(true);
243                                my_arx->set_statistics(1, V0, V0.rows()+1); 
244                        }
245                        else
246                        {
247                               
248                                V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
249                                V0(0,0) = 0;
250                                //V0(0,1) = -0.01;
251                                //V0(1,0) = -0.01;
252                                my_arx->set_constant(false);
253                                my_arx->set_statistics(1, V0, V0.rows()+1); 
254                               
255                        }                                               
256                       
257                        my_arx->set_parameters(window_size);
258                        my_arx->validate();
259
260                        previous_lognc = my_arx->posterior().lognc();
261
262                        /*
263                        vec mean = my_arx->posterior().mean();
264                        cout << mean << endl;
265                        */
266                }
267        }
268
269        ~model()
270        {
271                if(my_rarx!=NULL)
272                {
273                        delete my_rarx;
274                }
275                else
276                {
277                        delete my_arx;
278                }
279        }
280
281        void data_update(int time) 
282        {
283                vec data_vector;
284                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
285                { 
286                        data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
287                }
288                if(my_rarx!=NULL)
289                {       
290                        data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
291                        my_rarx->bayes(data_vector);
292                }
293                else
294                {
295                        vec pred_vec;
296                        pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
297                        my_arx->bayes(pred_vec,data_vector);
298                }
299        }
300
301        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to
302        {
303                vec condition_vector;
304                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
305                {
306                        condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1));
307                }
308
309                if(my_rarx!=NULL)
310                {
311                        pair<vec,mat> imp_samples = my_rarx->posterior->sample(sample_size,true);
312
313                        //cout << imp_samples.first << endl;                   
314                       
315                        vec sample_prediction;                 
316                        for(int t = 0;t<imp_samples.second.cols();t++)
317                        {
318                                vec lap_sample = condition_vector;
319                               
320                                if(has_constant)
321                                {
322                                        lap_sample.ins(lap_sample.size(),1.0);
323                                }
324                               
325                                lap_sample.ins(lap_sample.size(),(*LapRNG)());
326
327                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                             
328                        }
329
330                        return pair<vec,vec>(imp_samples.first,sample_prediction);
331                }       
332                else
333                {
334                        mat samples = my_arx->posterior().sample_mat(sample_size);                     
335                       
336                        vec sample_prediction;
337                        for(int t = 0;t<sample_size;t++)
338                        {
339                                vec gau_sample = condition_vector;
340                               
341                                if(has_constant)
342                                {
343                                        gau_sample.ins(gau_sample.size(),1.0);
344                                }
345                               
346                                gau_sample.ins(gau_sample.size(),randn());
347
348                                vec param_sample = samples.get_col(t);
349                                param_sample.set(param_sample.size()-1,sqrt(param_sample[param_sample.size()-1]));
350                                sample_prediction.ins(0,gau_sample*param_sample);                               
351                        }
352
353                        return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction);
354                }
355       
356        }
357
358
359        static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels)
360        {
361                set<set<pair<int,int>>> created_model_types;           
362
363                if(max_order == 1)//ukoncovacia vetva
364                {                       
365                        for(int channel = 0;channel<number_of_channels;channel++)//pre AR 1 model vytvori kombinace kanalov v prvom kroku poyadu
366                        {
367                                set<pair<int,int>> returned_type;
368                                returned_type.insert(pair<int,int>(channel,1));  //??
369                                created_model_types.insert(returned_type);
370                        }
371
372                        return created_model_types;
373                }
374                else
375                {
376                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);//tu uz mame ulozene kombinace o jeden krok dozadu  //rekuryivne volanie
377                        set<set<pair<int,int>>> returned_types;
378                       
379                        for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++)
380                        {                               
381                               
382                                for(int order = 1; order<=max_order; order++)
383                                {
384                                        for(int channel = 0;channel<number_of_channels;channel++)
385                                        {
386                                                set<pair<int,int>> returned_type;
387                                                pair<int,int> new_pair = pair<int,int>(channel,order);//??
388                                                if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) //??
389                                                {
390                                                        returned_type.insert((*model_ref).begin(),(*model_ref).end()); //co vlozi na zaciatok retuned_type?
391                                                        returned_type.insert(new_pair);                                                 
392
393                                                        returned_types.insert(returned_type);                                                   
394                                                }
395                                        }
396                                }
397                        }
398
399                        created_model_types.insert(returned_types.begin(),returned_types.end());
400
401                        return created_model_types;
402                }               
403        }
404};
405
406// ****************************************************
407//                           MAIN MAIN MAIN MAIN MAIN
408// ****************************************************
409int main ( int argc, char* argv[] ) 
410{
411       
412        // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from
413    // 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
414    // can be compared to the classical setup.
415       
416        itpp::Laplace_RNG LapRNG = Laplace_RNG();
417
418        char* folder_string  =  "C:\\RobustExperiments\\"; // "C:\\dataADClosePercDiff"; //
419        char* data_folder    = "data\\";
420        char* results_folder = "results\\";     
421
422        char dfstring[150];
423        strcpy(dfstring,folder_string);
424        strcat(dfstring,data_folder);
425        strcat(dfstring,commodity);
426        vector<char*> files = listFiles(dfstring);
427       
428        for(int contract=0;contract<files.size();contract++)
429        {
430                char *cdf_str = new char[strlen(dfstring) + 1];
431                strcpy(cdf_str, dfstring);             
432                strcat(cdf_str,files[contract]);       
433       
434                mat data_matrix;
435                ifstream myfile(cdf_str);
436                if (myfile.is_open())
437                {               
438                        string line;
439                        while(getline(myfile,line))
440                        {                       
441                                vec data_vector;
442                                while(line.find(',') != string::npos) //zmenil som ciarku za medzeru
443                                {
444                                        //line.erase(0,1); // toto som sem pridal
445                                        int loc2 = line.find('\n');
446                                        int loc  = line.find(',');
447                                        data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                           
448                                        line.erase(0,loc+1);                                   
449                                }
450
451                                data_matrix.ins_row(data_matrix.rows(),data_vector);
452                        }               
453
454                        myfile.close(); 
455                }
456                else
457                {
458                        cout << "Can't open data file!" << endl;
459                }
460       
461                //konec nacitavania dat
462                set<set<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,data_matrix.rows()); //volanie funkce kde robi kombinace modelov
463                                                                                                        //to priradime do model_types, data_matrix.row urcuje pocet kanalov dat
464                vector<model*> models;
465                for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++)
466                {// prechadza rozne typy kanalov, a poctu regresorov
467                        for(int window_size = max_window_size-1;window_size < max_window_size;window_size++)
468                        {
469                                models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy
470                                models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));
471                                models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix));
472                                models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));             
473                        }
474
475                        //set<pair<int,int>> empty_list;
476                        //models.push_back(new model(empty_list,false,true,100,0,&data_matrix));
477                }
478       
479                mat result_lognc;
480                // mat result_preds;
481
482                ofstream myfilew;
483                char rfstring[150];                                     
484                strcpy(rfstring,folder_string);
485                strcat(rfstring,results_folder);
486                strcat(rfstring,commodity);
487                strcat(rfstring,files[contract]);
488
489                /*
490                char predstring[150];                                   
491                strcpy(predstring,folder_string);
492                strcat(predstring,results_folder);
493                strcat(predstring,commodity);
494                strcat(predstring,"PRED");
495                strcat(predstring,files[contract]);
496                */
497
498                for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols()
499                {       
500                        cout << "Steps: " << time << endl;
501               
502                        if(time == max_model_order)
503                        { 
504                               
505                                myfilew.open(rfstring,ios::app);
506                                for(int i = 0;i<models.size();i++)
507                                {
508                                        for(set<pair<int,int>>::iterator ar_ref = models[i]->ar_components.begin();ar_ref != models[i]->ar_components.end();ar_ref++)
509                                        {
510                                                myfilew << (*ar_ref).second << (*ar_ref).first;                                                 
511                                        }
512
513                                        myfilew << ".";
514
515                                        if(models[i]->my_arx == NULL)
516                                        {
517                                                myfilew << "31";
518                                        }
519                                        else
520                                        {
521                                                myfilew << "30";
522                                        }
523                                       
524                                        if(models[i]->has_constant)
525                                        {
526                                                myfilew << "61";
527                                        }
528                                        else
529                                        {
530                                                myfilew << "60";
531                                        }
532
533                                        myfilew << ",999,999,999,";
534                                }
535
536                                myfilew << "888" << endl;
537                                myfilew.close();
538                        }
539               
540                        //pocet stlpcov data_matrix je pocet casovych krokov
541                        double cur_loglikelihood;
542                        // vec preds;
543                        int  prev_samples_nr;
544                        bool previous_switch = true;
545                        for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++)
546                        {//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
547                                (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne
548
549                                //cout << "Updated." << endl;                           
550
551                                if((*model_ref)->my_rarx!=NULL) //vklada normalizacni faktor do cur_res_lognc
552                                {
553                                        //cout << "Maxlik vertex(robust):" << (*model_ref)->my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
554                                        cur_loglikelihood = (*model_ref)->my_rarx->posterior->_ll();                           
555                                }
556                                else
557                                {
558                                        //cout << "Mean(classical):" << (*model_ref)->my_arx->posterior().mean() << endl;
559                                        double cur_lognc = (*model_ref)->my_arx->posterior().lognc();
560                                        cur_loglikelihood = cur_lognc-(*model_ref)->previous_lognc;
561                                                               
562                                        (*model_ref)->previous_lognc = cur_lognc;                               
563                                }
564                               
565                                /*
566                        if(time == max_window_size-1)
567                        {
568                                //***********************
569                                int sample_size = 100000;
570                                //***********************
571                               
572                                pair<vec,mat> samples;
573                                if((*model_ref)->my_arx!=NULL)
574                                {
575                                        mat samp_mat = (*model_ref)->my_arx->posterior().sample_mat(sample_size);
576                                        samples = pair<vec,mat>(ones(samp_mat.cols()),samp_mat);
577                                }
578                                else
579                                {
580                                        samples = (*model_ref)->my_rarx->posterior->sample(sample_size,true);                                   
581                                }
582
583                                char fstring[80];                                       
584                                strcpy(fstring,file_string);
585                                strcat(fstring,(*model_ref)->name);
586                                strcat(fstring,".txt");
587                               
588                                //cout << samples.first << endl;
589                               
590                                myfilew.open(fstring,ios::app);
591                               
592                               
593                                //for(int i = 0;i<samples.first.size();i++)
594                                //{
595                                //      myfilew << samples.first.get(i) << ",";
596                                //}
597                                //myfilew << endl;
598                               
599
600                                for(int j = 0;j<samples.second.rows()+1;j++)
601                                {                       
602                                        for(int i = 0;i<samples.second.cols();i++)
603                                        {
604                                                if(j!=samples.second.rows())
605                                                {
606                                                        myfilew << samples.second.get(j,i) << ",";
607                                                }
608                                               
609                                                //else
610                                                //{
611                                                //      myfilew << "0,";
612                                                //}
613                                               
614                                        }                               
615                                        myfilew << endl;
616                                }
617
618                                cout << "*************************************" << endl;
619
620                                myfilew.close();                               
621                        }
622                        */                             
623
624                                if(time>prediction_time)
625                                {                       
626                                        int samples_nr;                                 
627                                        if(previous_switch)
628                                        {
629                                                samples_nr = 10000;
630                                        }
631                                        else
632                                        {
633                                                samples_nr = prev_samples_nr;
634                                        }
635                                       
636                                       
637                                        // PREDICTIONS
638                                        pair<vec,vec> predictions = (*model_ref)->predict(samples_nr,time,&LapRNG);
639
640                                        /*
641                                        myfilew.open(predstring,ios::app);
642                                        for(int i=0;i<10000;i++)
643                                        {
644                                                if(i<predictions.second.size())
645                                                {
646                                                        myfilew << predictions.second.get(i) << ",";
647                                                }
648                                                else
649                                                {
650                                                        myfilew << ",";
651                                                }
652                                        }
653                                        myfilew << endl;
654                                        myfilew.close();
655                                        */
656
657                                        if(previous_switch)
658                                        {
659                                                prev_samples_nr = predictions.second.size();
660                                                samples_nr = prev_samples_nr;
661                                        }
662
663                                        previous_switch = !previous_switch;
664
665                                        double optimalInvestment = newtonRaphson(0,0.00001,predictions.second,utility_order);
666
667                                        if(abs(optimalInvestment)>max_investment)
668                                        {
669                                                optimalInvestment = max_investment*sign(optimalInvestment);
670                                        }
671
672
673                                        /*
674                                        vec utilityValues;
675                                        for(int j=0;j<1000;j++)
676                                        {
677                                                utilityValues.ins(utilityValues.length(),valueCRRAUtility(-0.5+0.001*j, predictions.second, utility_order));
678                                        }*/
679
680                                        double avg_prediction = (predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()));
681
682                                        (*model_ref)->predictions.ins((*model_ref)->predictions.size(),avg_prediction);
683                                               
684                                        myfilew.open(rfstring,ios::app);
685                               
686                                        /*
687                                        for(int j=0;j<utilityValues.length();j++)
688                                        {
689                                                myfilew << utilityValues.get(j) << ",";
690                                        }
691                                        myfilew << endl;
692                                        */
693
694                                        myfilew << avg_prediction << "," << optimalInvestment << "," << samples_nr << "," << cur_loglikelihood << ","; 
695                                        myfilew.close();
696                                }                       
697                        }
698
699                        if(time>prediction_time&&(time+1)<data_matrix.cols())
700                        {
701                                // REAL PRICE           
702                                myfilew.open(rfstring,ios::app);
703                                myfilew << data_matrix.get(0,time+1) << endl;
704                                myfilew.close();
705                        }       
706                }
707
708                for(vector<model*>::reverse_iterator model_ref = models.rbegin();model_ref!=models.rend();model_ref++)
709                {
710                        delete *model_ref;
711                }
712        }
713
714        return 0;
715}
716
717
Note: See TracBrowser for help on using the browser.