root/applications/robust/main.cpp @ 1423

/*!
\file
\brief Robust
\author Vasek Smidl 

 */

#include "estim/arx.h"
#include "robustlib.h"
#include <vector>
#include <iostream>
#include <fstream>
//#include <itpp/itsignal.h>
#include "windows.h"
#include "ddeml.h"
#include "stdio.h"
#include <dirent.h>
//#include <itpp/itoptim.h>

//#include "DDEClient.h"
//#include <conio.h>


using namespace itpp;
using namespace bdm;

//const int emlig_size = 2;
//const int utility_constant = 5;

const int max_model_order         = 2;
const double apriorno             = 0.01;
const int max_window_size         = 40;
const int utility_order           = 25;
const int prediction_time         = 30;
const double min_utility_argument = 0.001;
const double max_investment       = 3.0;
const int sample_size                     = 5000;
const char* commodity                 = "TY\\";

/*
HDDEDATA CALLBACK DdeCallback(
        UINT uType,     // Transaction type.
        UINT uFmt,      // Clipboard data format.
        HCONV hconv,    // Handle to the conversation.
        HSZ hsz1,       // Handle to a string.
        HSZ hsz2,       // Handle to a string.
        HDDEDATA hdata, // Handle to a global memory object.
        DWORD dwData1,  // Transaction-specific data.
        DWORD dwData2)  // Transaction-specific data.
{
        return 0;
}

void DDERequest(DWORD idInst, HCONV hConv, char* szItem)
{
        HSZ hszItem = DdeCreateStringHandle(idInst, szItem, 0);
        HDDEDATA hData = DdeClientTransaction(NULL,0,hConv,hszItem,CF_TEXT, 
                                                                        XTYP_ADVSTART,TIMEOUT_ASYNC , NULL); //TIMEOUT_ASYNC
        if (hData==NULL)
        {
                printf("Request failed: %s\n", szItem);
        }

        if (hData==0)
        {
                printf("Request failed: %s\n", szItem);
        }
}

DWORD WINAPI ThrdFunc( LPVOID n )
{    
    return 0;
}
*/
vector<char*> listFiles(char* dir){
        vector<char*> files;
                DIR *pDIR;
        struct dirent *entry;
        if( pDIR=opendir(dir) ){
                while(entry = readdir(pDIR)){
                        if( strcmp(entry->d_name, ".") != 0 && strcmp(entry->d_name, "..") != 0 )
                                                {                                                       
                                                        char *file_string = new char[strlen(entry->d_name) + 1];
                                                        strcpy(file_string, entry->d_name);                                                                                             
                                                        files.push_back(file_string);
                                                }
                }
                closedir(pDIR);
        }
                return files;
}

double valueCRRAUtility(const double &position, const pair<vec,vec> &samples, const int order)
{
        double value = 0;
        set<double> values;
                        
        for(int i=0;i<samples.first.length();i++)
        {
                double sample      = samples.second.get(i);
                double probability = samples.first.get(i);
                if((position*sample+1)>min_utility_argument)
                {
                        values.insert(probability*sample/pow(position*sample+1,order+1));
                }
                else
                {
                        values.insert(probability*(min_utility_argument-1)/position/pow(min_utility_argument,order+1));
                }
        }       
        
        for(set<double>::iterator val_ref = values.begin();val_ref!=values.end();val_ref++)
        {
                value+=(*val_ref);
        }

        return value;
}

double gradientCRRAUtility(const double &position, const pair<vec,vec> &samples, const int order)
{       
        double value = 0;
        set<double> values;

        for(int i=0;i<samples.first.length();i++)
        {
                double sample      = samples.second.get(i);
                double probability = samples.first.get(i);
                if((position*sample+1)>min_utility_argument)
                {
                        values.insert((-(order+1)*pow(sample,2)*probability)/pow(position*sample+1,order+2));
                }
        }

        for(set<double>::iterator val_ref = values.begin();val_ref!=values.end();val_ref++)
        {
                value+=(*val_ref);
        }

        return value;
}

double newtonRaphson(double startingPoint, double epsilon, pair<vec,vec> samples, int order)
{
        /*
        if(samples.length()>800)
        {
                samples.del(801,samples.size()-1);
        }
        */
        
        int count = 0;
        
        bool epsilon_reached = false;

        while(count<1000&&!epsilon_reached)
        {
                double cur_value = valueCRRAUtility(startingPoint,samples,order);
                double cur_gradient = gradientCRRAUtility(startingPoint,samples,order);

                startingPoint = startingPoint - cur_value/cur_gradient;

                if(cur_value<epsilon)
                {
                        epsilon_reached = true;
                }
        }

        if(count==100)
        {
                return startingPoint; // can be different!
        }
        else
        {
                return startingPoint;
        }
}

class model
{
public:
        set<pair<int,int>> ar_components;

        // Best thing would be to inherit the two models from a single souce, this is planned, but now structurally
        // problematic.
        RARX* my_rarx; //vzmenovane parametre pre triedu model
        ARXwin* my_arx;

        bool has_constant;
        int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit
        int  predicted_channel;
        mat* data_matrix;
        vec  predictions;
        char name[80];

        double previous_lognc;
        
        model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor
                  bool robust, 
                  bool has_constant, 
                  int window_size, 
                  int predicted_channel,
                  mat* data_matrix)
        {               
                this->ar_components.insert(ar_components.begin(),ar_components.end());

                strcpy(name,"M");

                for(set<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++)
                {
                        char buffer1[2];
                        char buffer2[2];
                        itoa((*ar_ref).first,buffer1,10);
                        itoa((*ar_ref).second,buffer2,10);

                        strcat(name,buffer1);
                        strcat(name,buffer2);
                        strcat(name,"_");
                }

                this->has_constant      = has_constant;

                if(has_constant)
                {
                        strcat(name,"C");
                }

                this->window_size       = window_size;
                this->predicted_channel = predicted_channel;
                this->data_matrix       = data_matrix;

                if(robust)
                {
                        previous_lognc = 0;
                        strcat(name,"R");

                        if(has_constant)
                        {
                                my_rarx = new RARX(ar_components.size()+1,window_size,true,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+3);
                                my_arx  = NULL;
                        }
                else
                        {
                                my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+2);
                                my_arx  = NULL;
                        }
                }
                else
                {
                        my_rarx = NULL;
                        my_arx  = new ARXwin();
                        mat V0;

                        if(has_constant)
                        {                               
                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst
                                V0(0,0) = 0;
                                my_arx->set_constant(true);
                                my_arx->set_statistics(1, V0, V0.rows()+1); 
                        }
                        else
                        {
                                
                                V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
                                V0(0,0) = 0;
                                //V0(0,1) = -0.01;
                                //V0(1,0) = -0.01;
                                my_arx->set_constant(false);
                                my_arx->set_statistics(1, V0, V0.rows()+1); 
                                
                        }                                               
                        
                        my_arx->set_parameters(window_size);
                        my_arx->validate();

                        previous_lognc = my_arx->posterior().lognc();

                        /*
                        vec mean = my_arx->posterior().mean();
                        cout << mean << endl;
                        */
                }
        }

        ~model()
        {
                if(my_rarx!=NULL)
                {
                        delete my_rarx;
                }
                else
                {
                        delete my_arx;
                }
        }

        void data_update(int time) 
        {
                vec data_vector;
                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
                {  
                        data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
                }
                if(my_rarx!=NULL)
                {       
                        data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
                        my_rarx->bayes(data_vector);
                }
                else
                {
                        vec pred_vec;
                        pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
                        my_arx->bayes(pred_vec,data_vector);
                }
        }

        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to
        {
                vec condition_vector;
                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
                {
                        condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1));
                }

                if(my_rarx!=NULL)
                {
                        pair<vec,mat> imp_samples = my_rarx->posterior->sample(sample_size,false);

                        //cout << imp_samples.first << endl;                    
                        
                        vec sample_prediction;                  
                        for(int t = 0;t<imp_samples.second.cols();t++)
                        {
                                vec lap_sample = condition_vector;
                                
                                if(has_constant)
                                {
                                        lap_sample.ins(lap_sample.size(),1.0);
                                }
                                
                                lap_sample.ins(lap_sample.size(),(*LapRNG)());

                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                              
                        }

                        return pair<vec,vec>(imp_samples.first,sample_prediction);
                }       
                else
                {
                        mat samples = my_arx->posterior().sample_mat(sample_size);                      
                        
                        vec sample_prediction;
                        for(int t = 0;t<sample_size;t++)
                        {
                                vec gau_sample = condition_vector;
                                
                                if(has_constant)
                                {
                                        gau_sample.ins(gau_sample.size(),1.0);
                                }
                                
                                gau_sample.ins(gau_sample.size(),randn());

                                vec param_sample = samples.get_col(t);
                                param_sample.set(param_sample.size()-1,sqrt(param_sample[param_sample.size()-1]));
                                sample_prediction.ins(0,gau_sample*param_sample);                               
                        }

                        return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction);
                }
        
        }


        static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels)
        {
                set<set<pair<int,int>>> created_model_types;            

                if(max_order == 1)//ukoncovacia vetva
                {                       
                        for(int channel = 0;channel<number_of_channels;channel++)//pre AR 1 model vytvori kombinace kanalov v prvom kroku poyadu
                        {
                                set<pair<int,int>> returned_type;
                                returned_type.insert(pair<int,int>(channel,1));  //??
                                created_model_types.insert(returned_type);
                        }

                        return created_model_types;
                }
                else
                {
                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);//tu uz mame ulozene kombinace o jeden krok dozadu  //rekuryivne volanie
                        set<set<pair<int,int>>> returned_types;
                        
                        for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++)
                        {                               
                                
                                for(int order = 1; order<=max_order; order++)
                                {
                                        for(int channel = 0;channel<number_of_channels;channel++)
                                        {
                                                set<pair<int,int>> returned_type;
                                                pair<int,int> new_pair = pair<int,int>(channel,order);//??
                                                if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) //??
                                                {
                                                        returned_type.insert((*model_ref).begin(),(*model_ref).end()); //co vlozi na zaciatok retuned_type?
                                                        returned_type.insert(new_pair);                                                 

                                                        returned_types.insert(returned_type);                                                   
                                                }
                                        }
                                }
                        }

                        created_model_types.insert(returned_types.begin(),returned_types.end());

                        return created_model_types;
                }               
        }
};

// ****************************************************
//                           MAIN MAIN MAIN MAIN MAIN
// ****************************************************
int main ( int argc, char* argv[] ) 
{
        
        // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from 
    // 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
    // can be compared to the classical setup.
        
        itpp::Laplace_RNG LapRNG = Laplace_RNG();

        char* folder_string  =  "C:\\RobustExperiments\\"; // "C:\\dataADClosePercDiff"; //
        char* data_folder    = "data\\";
        char* results_folder = "results\\";     

        char dfstring[150];
        strcpy(dfstring,folder_string);
        strcat(dfstring,data_folder);
        strcat(dfstring,commodity);
        vector<char*> files = listFiles(dfstring);
        
        for(int contract=0;contract<files.size();contract++)
        {
                char *cdf_str = new char[strlen(dfstring) + 1];
                strcpy(cdf_str, dfstring);              
                strcat(cdf_str,files[contract]);        
        
                mat data_matrix;
                ifstream myfile(cdf_str);
                if (myfile.is_open())
                {               
                        string line;
                        while(getline(myfile,line))
                        {                       
                                vec data_vector;
                                while(line.find(',') != string::npos) //zmenil som ciarku za medzeru
                                {
                                        //line.erase(0,1); // toto som sem pridal
                                        int loc2 = line.find('\n');
                                        int loc  = line.find(',');
                                        data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                           
                                        line.erase(0,loc+1);                                    
                                }

                                data_matrix.ins_row(data_matrix.rows(),data_vector);
                        }               

                        myfile.close(); 
                }
                else
                {
                        cout << "Can't open data file!" << endl;
                }
        
                //konec nacitavania dat
                set<set<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,data_matrix.rows()); //volanie funkce kde robi kombinace modelov
                                                                                                        //to priradime do model_types, data_matrix.row urcuje pocet kanalov dat
                vector<model*> models;
                for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++)
                {// prechadza rozne typy kanalov, a poctu regresorov
                        for(int window_size = max_window_size-1;window_size < max_window_size;window_size++)
                        {
                                //if(model_type->size()<max_model_order)
                                //{
                                        models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy
                                        models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));
                                //}
                                models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix));
                                models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));              
                        }

                        //set<pair<int,int>> empty_list;
                        //models.push_back(new model(empty_list,false,true,100,0,&data_matrix));
                }
        
                mat result_lognc;
                // mat result_preds;

                ofstream myfilew;
                char rfstring[150];                                     
                strcpy(rfstring,folder_string);
                strcat(rfstring,results_folder);
                strcat(rfstring,commodity);
                strcat(rfstring,files[contract]);

                /*
                char predstring[150];                                   
                strcpy(predstring,folder_string);
                strcat(predstring,results_folder);
                strcat(predstring,commodity);
                strcat(predstring,"PRED");
                strcat(predstring,files[contract]);
                */

                for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols() 
                {       
                        cout << "Steps: " << time << endl;
                
                        if(time == max_model_order)
                        { 
                                /*
                                myfilew.open(rfstring,ios::app);
                                for(int i = 0;models.size();i++)
                                {
                                        for(set<pair<int,int>>::iterator ar_ref = models[i]->ar_components.begin();ar_ref != models[i]->ar_components.end();ar_ref++)
                                        {
                                                myfilew << (*ar_ref).second << (*ar_ref).first;                                                 
                                        }

                                        myfilew << ".";

                                        if(models[i]->my_arx == NULL)
                                        {
                                                myfilew << "31";
                                        }
                                        else
                                        {
                                                myfilew << "30";
                                        }
                                        
                                        if(models[i]->has_constant)
                                        {
                                                myfilew << "61";
                                        }
                                        else
                                        {
                                                myfilew << "60";
                                        }

                                        myfilew << ",999,999,999,";
                                }

                                myfilew << "888" << endl;
                                myfilew.close();
                                */
                        }
                
                        //pocet stlpcov data_matrix je pocet casovych krokov
                        double cur_loglikelihood;
                        // vec preds;
                        int  prev_samples_nr;
                        bool previous_switch = true;
                        for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++) //.begin()+1;model_ref++)
                        {//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
                                (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne

                                //cout << "Updated." << endl;                           

                                if((*model_ref)->my_rarx!=NULL) //vklada normalizacni faktor do cur_res_lognc
                                {
                                        //cout << "Maxlik vertex(robust):" << (*model_ref)->my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
                                        cur_loglikelihood = (*model_ref)->my_rarx->posterior->_ll();                            
                                }
                                else
                                {
                                        //cout << "Mean(classical):" << (*model_ref)->my_arx->posterior().mean() << endl;
                                        double cur_lognc = (*model_ref)->my_arx->posterior().lognc();
                                        cur_loglikelihood = cur_lognc-(*model_ref)->previous_lognc;
                                                                
                                        (*model_ref)->previous_lognc = cur_lognc;                               
                                }
                                
                                /*
                        if(time == max_window_size-1)
                        {
                                //***********************
                                int sample_size = 100000;
                                //***********************
                                
                                pair<vec,mat> samples;
                                if((*model_ref)->my_arx!=NULL)
                                {
                                        mat samp_mat = (*model_ref)->my_arx->posterior().sample_mat(sample_size);
                                        samples = pair<vec,mat>(ones(samp_mat.cols()),samp_mat);
                                }
                                else
                                {
                                        samples = (*model_ref)->my_rarx->posterior->sample(sample_size,true);                                   
                                }

                                char fstring[80];                                       
                                strcpy(fstring,file_string);
                                strcat(fstring,(*model_ref)->name);
                                strcat(fstring,".txt");
                                
                                //cout << samples.first << endl;
                                
                                myfilew.open(fstring,ios::app);
                                
                                
                                //for(int i = 0;i<samples.first.size();i++)
                                //{
                                //      myfilew << samples.first.get(i) << ",";
                                //}
                                //myfilew << endl;
                                

                                for(int j = 0;j<samples.second.rows()+1;j++)
                                {                       
                                        for(int i = 0;i<samples.second.cols();i++)
                                        {
                                                if(j!=samples.second.rows())
                                                {
                                                        myfilew << samples.second.get(j,i) << ",";
                                                }
                                                
                                                //else
                                                //{
                                                //      myfilew << "0,";
                                                //}
                                                
                                        }                               
                                        myfilew << endl;
                                }

                                cout << "*************************************" << endl;

                                myfilew.close();                                
                        }
                        */                              

                                if(time>prediction_time)
                                {                       
                                        int samples_nr;                                 
                                        if(previous_switch)
                                        {
                                                samples_nr = sample_size;
                                        }
                                        else
                                        {
                                                samples_nr = prev_samples_nr;
                                        }
                                        
                                        
                                        // PREDICTIONS
                                        pair<vec,vec> predictions = (*model_ref)->predict(samples_nr,time,&LapRNG);

                                        /*
                                        myfilew.open(predstring,ios::app);
                                        for(int i=0;i<10000;i++)
                                        {
                                                if(i<predictions.second.size())
                                                {
                                                        myfilew << predictions.second.get(i) << ",";
                                                }
                                                else
                                                {
                                                        myfilew << ",";
                                                }
                                        }
                                        myfilew << endl;
                                        myfilew.close();
                                        */

                                        if(previous_switch)
                                        {
                                                prev_samples_nr = predictions.second.size();
                                                samples_nr = prev_samples_nr;
                                        }

                                        previous_switch = !previous_switch;

                                        double optimalInvestment = newtonRaphson(0,0.00001,predictions,utility_order);

                                        if(abs(optimalInvestment)>max_investment)
                                        {
                                                optimalInvestment = max_investment*sign(optimalInvestment);
                                        }


                                        /*
                                        vec utilityValues;
                                        for(int j=0;j<1000;j++)
                                        {
                                                utilityValues.ins(utilityValues.length(),valueCRRAUtility(-0.5+0.001*j, predictions.second, utility_order));
                                        }*/

                                        double avg_prediction = (predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()));

                                        (*model_ref)->predictions.ins((*model_ref)->predictions.size(),avg_prediction);
                                                
                                        myfilew.open(rfstring,ios::app);
                                
                                        /*
                                        for(int j=0;j<utilityValues.length();j++)
                                        {
                                                myfilew << utilityValues.get(j) << ",";
                                        }
                                        myfilew << endl;
                                        */

                                        myfilew << avg_prediction << "," << optimalInvestment << "," << samples_nr << "," << cur_loglikelihood << ","; 
                                        myfilew.close();
                                }                       
                        }

                        if(time>prediction_time&&(time+1)<data_matrix.cols())
                        {
                                // REAL PRICE           
                                myfilew.open(rfstring,ios::app);
                                myfilew << data_matrix.get(0,time+1) << endl;
                                myfilew.close();
                        }       
                }

                for(vector<model*>::reverse_iterator model_ref = models.rbegin();model_ref!=models.rend();model_ref++)
                {
                        delete *model_ref;
                }
        }

        return 0;
}