Context Navigation

main.cpp @ 1412

Revision 1412, 18.3 kB (checked in by sindj, 13 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
25	using namespace itpp;
26	using namespace bdm;
27
28	//const int emlig_size = 2;
29	//const int utility_constant = 5;
30
31	const int max_model_order = 2;
32	const double apriorno = 0.01;
33	const int max_window_size = 40;
34	const int utility_order = 25;
35	const int prediction_time = 30;
36	const double min_utility_argument = 0.001;
37	const double max_investment = 10.0;
38	const char* commodity = "CL\\";
39
40	/*
41	HDDEDATA 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
54	void 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
70	DWORD WINAPI ThrdFunc( LPVOID n )
71	{
72	return 0;
73	}
74	*/
75	vector<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
93	double 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
113	double 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(positionsample+1,order+2);
123	}
124	}
125
126	return value;
127	}
128
129	double 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
163	class model
164	{
165	public:
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(2apriorno),sqrt(2apriorno),ar_components.size()+3);
224	my_arx = NULL;
225	}
226	else
227	{
228	my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2apriorno),sqrt(2apriorno),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	// ****************************************************
409	int 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.95y_(t-1)+0.05y_(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.firstpredictions.second)/(predictions.firstones(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.

Context Navigation

root/applications/robust/main.cpp @ 1412

Download in other formats: