Index: applications/robust/main.cpp
===================================================================
--- applications/robust/main.cpp (revision 1406)
+++ applications/robust/main.cpp (revision 1412)
@@ -16,5 +16,6 @@
 #include "ddeml.h"
 #include "stdio.h"
-#include <itpp/itoptim.h>
+#include <dirent.h>
+//#include <itpp/itoptim.h>
 
 //#include "DDEClient.h"
@@ -28,9 +29,12 @@
 //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   = 19;
-const int prediction_time = 30;
+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       = 10.0;
+const char* commodity		      = "CL\\";
 
 /*
@@ -69,4 +73,21 @@
 }
 */
+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 vec &samples, const int order)
@@ -77,5 +98,12 @@
 	{
 		double sample = samples.get(i);
-		value += sample/pow(position*sample+1,order+1);
+		if((position*sample+1)>min_utility_argument)
+		{
+			value += sample/pow(position*sample+1,order+1);
+		}
+		else
+		{
+			value += (min_utility_argument-1)/position/pow(min_utility_argument,order+1);
+		}
 	}	
 	
@@ -90,5 +118,8 @@
 	{
 		double sample = samples.get(i);
-		value += (-(order+1)*pow(sample,2))/pow(position*sample+1,order+2);
+		if((position*sample+1)>min_utility_argument)
+		{
+			value += (-(order+1)*pow(sample,2))/pow(position*sample+1,order+2);
+		}
 	}
 
@@ -129,10 +160,4 @@
 	}
 }
-	
-	
-
-	
-
-
 
 class model
@@ -242,4 +267,16 @@
 	}
 
+	~model()
+	{
+		if(my_rarx!=NULL)
+		{
+			delete my_rarx;
+		}
+		else
+		{
+			delete my_arx;
+		}
+	}
+
 	void data_update(int time) 
 	{
@@ -309,5 +346,7 @@
 				gau_sample.ins(gau_sample.size(),randn());
 
-				sample_prediction.ins(0,gau_sample*samples.get_col(t));				
+				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);				
 			}
 
@@ -365,7 +404,7 @@
 };
 
-
-
-
+// ****************************************************
+//			     MAIN MAIN MAIN MAIN MAIN
+// ****************************************************
 int main ( int argc, char* argv[] ) 
 {
@@ -375,114 +414,154 @@
     // can be compared to the classical setup.
 	
-	itpp::Laplace_RNG LapRNG = Laplace_RNG();	
-	
-	vector<vector<string>> strings;
-
-	char* file_string =  "C:\\CD2"; // "C:\\dataADClosePercDiff"; //  
-
-	char dfstring[80];
-	strcpy(dfstring,file_string);
-	strcat(dfstring,".txt");
-	
-	
-	mat data_matrix;
-	ifstream myfile(dfstring);
-	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++)
-		{
-			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 fstring[80];					
-	strcpy(fstring,file_string);
-	//strcat(fstring,"lognc.txt");
-	strcat(fstring,"preds.txt");	
-
-	for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols() 
-	{ 	
-		cout << "Steps: " << time << endl;
+	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++)
+			{
+				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_window_size)
-		{
-			exit(1);
-		}*/
+			if(time == max_model_order)
+			{ 
+				
+				myfilew.open(rfstring,ios::app);
+				for(int i = 0;i<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
-		vec cur_res_lognc;
-		// vec preds;
-		vector<string> nazvy;
-		for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();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 (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu?
-
-			if((*model_ref)->my_rarx!=NULL)	//vklada normalizacni faktor do cur_res_lognc
-			{
-				//cout << "Maxlik vertex:" << (*model_ref)->my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
-				cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_rarx->posterior->_ll());				
-			}
-			else
-			{
-				double cur_lognc = (*model_ref)->my_arx->posterior().lognc();
-				double cur_ll = cur_lognc-(*model_ref)->previous_lognc;
+			//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++)
+			{//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<(*model_ref)->window_size)
-				{
-					cur_ll-=1.83787706640935;
-				}*/
-				
-				(*model_ref)->my_arx->_ll();
-
-				cur_res_lognc.ins(cur_res_lognc.size(),cur_ll);
-				
-				(*model_ref)->previous_lognc = cur_lognc;				
-			}
-			
-			/*
 			if(time == max_window_size-1)
 			{
@@ -541,278 +620,95 @@
 				myfilew.close();				
 			}
-			*/		
-
-			if(time>prediction_time)
-			{			
-				// PREDICTIONS
-				pair<vec,vec> predictions = (*model_ref)->predict(5000,time,&LapRNG);
-
-				/*
-				cout << predictions.first << endl << endl << predictions.second << endl << "*************************************" ;
-				pause(5);
-				*/	
-
-				double optimalInvestment = newtonRaphson(0,0.00001,predictions.second,utility_order);
-
-				/*
-				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);
+			*/				
+
+				if(time>prediction_time)
+				{			
+					int samples_nr;					
+					if(previous_switch)
+					{
+						samples_nr = 10000;
+					}
+					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.second,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(fstring,ios::app);
-				
-				/*
-				for(int j=0;j<utilityValues.length();j++)
-				{
-					myfilew << utilityValues.get(j) << ",";
-				}
-				myfilew << endl;
-				*/
-
-				myfilew << avg_prediction << "," << optimalInvestment << ","; 
+					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();
-			}
-			
-
-			//preds.ins(0,data_matrix.get(0,time+1));
-		}
-
-		if(time>prediction_time)
-		{
-			// REAL PRICE		
-			myfilew.open(fstring,ios::app);
-			myfilew << data_matrix.get(0,time+1) << endl;
-			myfilew.close();
-		}
-		
-		
-
-		result_lognc.ins_col(result_lognc.cols(),cur_res_lognc);
-		//result_preds.ins_col(result_preds.cols(),preds);
-
-		// cout << "Updated." << endl;	
-							
-		/*
-		myfilew.open(fstring,ios::app);
-		
-		// myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
-		
-		if(time == max_model_order)
-		{ 
-			for(int i = 0;i<cur_res_lognc.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 << "1";
-				}
-				else
-				{
-					myfilew << "0";
-				}
-					
-				if(models[i]->has_constant)
-				{
-					myfilew << "1";
-				}
-				else
-				{
-					myfilew << "0";
-				}
-
-				myfilew << ",";
-			}
-
-			myfilew << endl;
-		}
-		
-		
-		for(int i = 0;i<cur_res_lognc.size();i++)
-		{
-		   myfilew << cur_res_lognc[i] << ' ';//zmenil som ciarku ze medzeru
-		}
-
-		myfilew << endl;		
-		myfilew.close();	
-		*/
-}
-
-	// EXPERIMENT: One step ahead price prediction. Comparison of classical and robust model using optimal trading
-    //             with maximization of logarithm of one-step ahead wealth.
-
-	
-		
-		/*
-		cout << "One experiment finished." << endl;
-
-		ofstream myfile;
-		myfile.open("c:\\robust_ar1.txt",ios::app);
-		myfile << endl;
-		myfile.close();
-
-		myfile.open("c:\\robust_ar2.txt",ios::app);
-		myfile << endl;
-		myfile.close();*/
-	
-
-	//emlig* emlig1 = new emlig(emlig_size);
-
-	//emlig1->step_me(0);
-	//emlig* emlig2 = new emlig(emlig_size);
-	
-	/*
-	emlig1->set_correction_factors(4);
-
-	for(int j = 0;j<emlig1->correction_factors.size();j++)
-	{
-		for(set<my_ivec>::iterator vec_ref = emlig1->correction_factors[j].begin();vec_ref!=emlig1->correction_factors[j].end();vec_ref++)
-		{
-			cout << j << "    ";
-			
-			for(int i=0;i<(*vec_ref).size();i++)
-			{
-				cout << (*vec_ref)[i];
-			}
-
-			cout << endl;
-		}
-	}*/
-	
-	/*
-    vec condition5 = "1.0 1.0 1.01";//"-0.3 1.7 1.5";
-
-	emlig1->add_condition(condition5);
-	//emlig1->step_me(0);
-
-
-	vec condition1a = "-1.0 1.02 0.5";
-	//vec condition1b = "1.0 1.0 1.01";
-	emlig1->add_condition(condition1a);
-	//emlig2->add_condition(condition1b);
-
-	vec condition2a = "-0.3 1.7 1.5";
-	//vec condition2b = "-1.0 1.0 1.0";
-	emlig1->add_condition(condition2a);
-	//emlig2->add_condition(condition2b);
-
-	vec condition3a = "0.5 -1.01 1.0";
-	//vec condition3b = "0.5 -1.01 1.0";
-
-	emlig1->add_condition(condition3a);
-	//emlig2->add_condition(condition3b);	
-
-	vec condition4a = "-0.5 -1.0 1.0";
-	//vec condition4b = "-0.5 -1.0 1.0";	
-
-	emlig1->add_condition(condition4a);
-	//cout << "************************************************" << endl;
-	//emlig2->add_condition(condition4b);
-	//cout << "************************************************" << endl;
-	
-	//cout << emlig1->minimal_vertex->get_coordinates();
-	
-	//emlig1->remove_condition(condition3a);
-	//emlig1->step_me(0);
-	//emlig1->remove_condition(condition2a);
-	//emlig1->remove_condition(condition1a);
-	//emlig1->remove_condition(condition5);
-	
-
-	//emlig1->step_me(0);
-	//emlig2->step_me(0);
-	
-
-	// DA SE POUZIT PRO VYPIS DO SOUBORU
-	// emlig1->step_me(0);
-
-	//emlig1->remove_condition(condition1);
-	
-	
-
-	
-	
-	/*
-	for(int i = 0;i<100;i++)
-	{
-		cout << endl << "Step:" << i << endl;		
-
-		double condition[emlig_size+1];		
-
-		for(int k = 0;k<=emlig_size;k++)
-		{
-			condition[k] = (rand()-RAND_MAX/2)/1000.0;		
-		}
-			
-
-		vec* condition_vec = new vec(condition,emlig_size+1);
-		emlig1->add_condition(*condition_vec);
-
-		/*
-		for(polyhedron* toprow_ref = emlig1->statistic.rows[emlig_size]; toprow_ref != emlig1->statistic.end_poly; toprow_ref = toprow_ref->next_poly)
-		{
-			cout << ((toprow*)toprow_ref)->probability << endl;
-		}
-		*/
-		/*
-		cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
-	
-		/*
-		if(i-emlig1->number_of_parameters >= 0)
-		{
-			pause(30);
-		}
-		*/
-
-		// emlig1->step_me(i);
-		
-		/*
-		vector<int> sizevector;
-		for(int s = 0;s<=emlig1->number_of_parameters;s++)
-		{
-			sizevector.push_back(emlig1->statistic_rowsize(s));
-		}
-		*/
-	//}
-    
-
-
-	
-	/*
-	emlig1->step_me(1);
-
-	vec condition = "2.0 0.0 1.0";	
-
-	emlig1->add_condition(condition);
-
-	vector<int> sizevector;
-	for(int s = 0;s<=emlig1->number_of_parameters;s++)
-	{
-		sizevector.push_back(emlig1->statistic_rowsize(s));
-	}
-
-	emlig1->step_me(2);
-
-	condition = "2.0 1.0 0.0";
-
-	emlig1->add_condition(condition);
-
-	sizevector.clear();
-	for(int s = 0;s<=emlig1->number_of_parameters;s++)
-	{
-		sizevector.push_back(emlig1->statistic_rowsize(s));
-	}
-	*/
+			}	
+		}
+
+		for(vector<model*>::reverse_iterator model_ref = models.rbegin();model_ref!=models.rend();model_ref++)
+		{
+			delete *model_ref;
+		}
+	}
 
 	return 0;
