Context Navigation

main.cpp @ 1423

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

Note: See TracBrowser for help on using the browser.

Context Navigation

root/applications/robust/main.cpp @ 1423

Download in other formats: