FE587
Home
Description
Programs

FE 587 - C/C++ Programming for Financial Engineering

Final on 26.12.2015

Assignment

  • Implement the Particle filter starting from the following skeleton

  • Final practice: Extend the Kalman filter implementation to a 2 x 2 transition matrix.

Programs covered on 26.12.2015

File Output
#include <iostream>
#include <fstream>
#include <ctime>
#include <cstdlib>


using namespace std;

int main() {
	 ofstream ofs;
	 //my_output_file.open ("test.txt", std::ofstream::out | std::ofstream::app);
	 ofs.open ("test.txt");

	 for (int i=0; i<100; i++) {
		 ofs << time(NULL) << " " << rand()%1000 << endl;
	 }

	 ofs.close();

	 return 0;
}
File Input
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

using namespace std;

int main() {
	 vector<int> time_stamps;
	 vector<double> data;

	 ifstream ifs;
	 double x;
	 int ts;
	 ifs.open ("../files/test.txt");

	 while (not ifs.eof()) {
		 ifs >> ts;
		 ifs >> x;
		 time_stamps.push_back(ts);
		 data.push_back(x);
	 }

	 ifs.close();

	 for (int i=0; i<data.size(); i++) {
		 cout << data[i] << endl;
	 }

	 return 0;

}

Programs covered on 19.12.2015

Particle Filter

#include <iostream>
#include <cmath>

#include <vector>

using namespace std;

class ParticleFilter {
	int N;  // Number of particles
	double A; // Transition Matrix
	double C; // Observation Matrix
	double Q; // Transition Covariance
	double R; // Observation Covariance
public:
	ParticleFilter(int _N = 100, double _A = 1, double _C = 1, double _Q = 1, double _R = 1) {
		A = _A;
		C = _C;
		Q = _Q;
		R = _R;
		N = _N;
	};

	void resample(vector<double>& W, vector<double>& x) {

	}

	double filter(vector<double>& y, double m, double Sig, vector<double>& m_f, vector<double>& Cov_f) {
		vector<double> x(N); // Particle coordinate
		vector<double> Weight(N);
		double e;
		for (int j=0; j<N; j++) {
			x[j] = m + sqrt(Sig)*rand_normal();
			Weight[j] = 1.0/N;
		}

		for (unsigned int t = 0; t<y.size(); t++) {
			for (int j=0; j<N; j++) {
				e = y[t] - C*x[j];
				Weight[j] = exp( -0.5*log(2*M_PI*R) - 0.5*e*e/R);
			}

			resample(Weight, x);

			// Predict
			for (int j=0; j<N; j++) {
				x[j] = A*x[j] + sqrt(Q)*rand_normal();
			}


		}

		return ll;

	};
};



int main() {
	ParticleFilter pf(500, 1, 0.01, 0.05);
	double yy[] = {0.3,0.7,0.95,3.2,5.0};
	vector<double> y(yy, yy+sizeof(yy)/sizeof(double));
	vector<double> m(y.size());
	vector<double> S(y.size());

	pf.filter(y, 0, 50, m, S);

	for (unsigned int i=0; i<y.size(); i++) {
		cout << y[i] << '\t' << m[i] << '\t' << S[i] << endl;
	}

	return 0;
}
Kalman Filter

#include <iostream>
#include <cmath>

#include <vector>

using namespace std;

class KalmanFilter {
	double A; // Transition Matrix
	double C; // Observation Matrix
	double Q; // Transition Covariance
	double R; // Observation Covariance
public:
	KalmanFilter(double _A = 1, double _C = 1, double _Q = 1, double _R = 1) {
		A = _A;
		C = _C;
		Q = _Q;
		R = _R;
	};

	double filter(vector<double>& y, double m, double Sig, vector<double>& m_f, vector<double>& Cov_f) {
		double Gain;
		double ll = 0;
		double e = 0, mu, S;

		for (unsigned int i = 0; i<y.size(); i++) {
			mu = C*m;
			S = C*Sig*C + R;
			e = (y[i]-mu);
			ll += -0.5*log(2*M_PI*S) - 0.5*e*e/S;

			// Update

			Gain = Sig*C/(C*Sig*C+R);
			m = m + Gain*(y[i] - C*m);
			Sig = Sig - Gain*C*Sig;

			m_f[i] = m;
			Cov_f[i] = Sig;

			// Predict
			m = A*m;
			Sig = A*Sig*A + Q;

		}

		return ll;

	};
};



int main() {
	KalmanFilter kf(500, 1, 0.01, 0.05);
	double yy[] = {0.3,0.7,0.95,3.2,5.0};
	vector<double> y(yy, yy+sizeof(yy)/sizeof(double));
	vector<double> m(y.size());
	vector<double> S(y.size());

	double ll = kf.filter(y, 0, 50, m, S);

	for (unsigned int i=0; i<y.size(); i++) {
		cout << y[i] << '\t' << m[i] << '\t' << S[i] << endl;
	}

	cout << "LogLikelihood = " << ll << endl;
	return 0;
}

Midterm on 12.12.2015

Programs covered on 05.12.2015

#include <iostream>
#include <vector>
using namespace std;


typedef vector<double> Row;


class feMatrix {
	int nRows;
	int nCols;
	vector<Row> data;

public:

	feMatrix(int nr = 1, int nc = 1, double def_val = 0.0) {
		nRows = nr;
		nCols = nc;
		data.resize(nRows);
		for (int i=0; i<nRows; i++) {
			data[i].resize(nCols);
		}

	}

	Row& operator[](int i) {
		return data[i];
	}

	int size(int dim = 1) {
		if (dim==1) return nRows;
		else return nCols;
	}
};

ostream& operator<<(ostream& os, Row& x) {
	for (int i=0; i<x.size(); i++) {
		os << x[i] << "\t";
	}
	return os;
}


ostream& operator<<(ostream& os, feMatrix& x) {
	for (int i=0; i<x.size(1); i++) {
		os << x[i] << endl;
	}
	return os;
}


void conditional(feMatrix& mu, feMatrix& Sigma, feMatrix& y, feMatrix& cond_mu, feMatrix& cond_V) {

}

int main() {
	feMatrix Sigma(2,2);
	feMatrix mu(2,1);
	feMatrix y(1,1);

	feMatrix cond_Sigma(1,1);
	feMatrix cond_mu(1,1);


	Sigma[0][0] = 1.0;
	Sigma[0][1] = 0.3;
	Sigma[1][0] = 0.3;
	Sigma[1][1] = 1.0;

	mu[0][0] = 0.0;
	mu[1][1] = 0.0;

	y[0][0] = 3.3;

	cout << Sigma;
	cout << mu;

	conditional(mu, Sigma, y, cond_mu, cond_Sigma);


}

Assignment due 5.12.2015

  1. (Harder) Extend the definition of the Vect object such that the following code segment works

#include <iostream>
#include <vector>

using namespace std;


int main ()
{
  // constructors used in the same order as described above:
  Vect<int> first;                                // empty Vect of ints
  Vect<int> second (4,100);                       // four ints with value 100
  Vect<int> third (second.begin(),second.end());  // iterating through second
  Vect<int> fourth (third);                       // a copy of third

  // the iterator constructor can also be used to construct from arrays:
  int myints[] = {16,2,77,29, 691, 72, 7892, 8943};
  //Vect&lt;<span class="builtin">int</span>&gt; fifth (myints+2, myints + sizeof(myints) / sizeof(<span class="builtin">int</span>) );
  Vect<int> fifth (myints+2, myints + 5 );

  cout << "The contents of fifth are:";
  for (Vect<int>::iterator it = fifth.begin(); it != fifth.end(); ++it)
    cout << ' ' << *it;
  cout << '\n';

  return 0;
}

  1. Make the following code work with a vector object instead of Vect – you can change method names 

int main() {
	Vect<complex_num> v(5, complex_num(-1,3));
	cout << v << endl;

	v.extend(30);
	cout << v << endl;
	return 0;
}

Programs covered on 28.11.2015

template example
#include <iostream>
using namespace std;


struct complex_num {
	double re;
	double im;
	complex_num(double r=0, double i=0) {
		re = r;
		im = i;
	}
};


ostream& operator<<(ostream& os, complex_num& c) {
	os << c.re << (c.im < 0 ? "" : "+") << c.im << "j";
	return os;
}

template <typename T>
class obj {
	T* data;
	int N;

public:
	obj(int L=10, T def_value=0) {
		data = new T[L];
		for (int i=0; i<L; i++) data[i] = def_value;
		N = L;
	}

	void print() {
		for (int i=0; i<N; i++) cout << data[i] << " ";
		cout << endl;
	}

	~obj() {
		if (data) delete data;
	}
};



int main() {

	obj<double> v(5, -3);
	obj<char> v2(5, 'X');
	obj<complex_num> v3(5, complex_num(-1,3));

	v.print();
	v2.print();
	v3.print();

	return 0;
}
Class Example : Vect object
#include <iostream>
using namespace std;

class Vect {
	double* data;
	int N;
	int Reserved;

public:
	Vect(int L=10) {
		if (L>0)
			Reserved = L*4;
			data = new double[Reserved];
			for (int i=0;i<L; i++) data[i] = 2*i;
			N = L;
	}

	~Vect() {
		if (data!=NULL) delete data;
	}

	void extend(int M = 10) {
		if ((N+M) < Reserved) {
			for (int i=0; i<M; i++) data[i+N] = -1;
			N+= M;
		}
		else {
			double* temp = new double[2*Reserved+M];
			if (temp!=NULL) {
				Reserved = 2*Reserved+M;
				for (int i=0; i<N; i++) temp[i] = data[i];
				for (int i=0; i<M; i++) temp[i+N] = -1;

				N += M;
				delete data;
				data = temp;
			}
		}
		return;
	}

	int Length() {
		return N;
	}

	double* Data() {
		return data;
	}
};

ostream& operator<<(ostream& os, Vect& v) {
	double *d = v.Data();
	for (int i=0; i<v.Length(); i++) {
		os << d[i] << " ";
	}
	return os;
}


int main() {
	Vect v(5);
	cout << v << endl;

	v.extend(30);
	cout << v << endl;
	return 0;
}

Optional Assignment for 28.11.2015: Extend the ts (time-series) object with a weighted average where the weights are provided as another ts object.

Time Series Object
#include <cstdlib>
#include <iostream>
#include <ctime>
using namespace std;

// Time series object
struct ts {
	int t_begin;  // time index of the first sample
	int t_end;    // time index of the last sample
	double* data; // Pointer to data buffer
	int N;        // Number of samples, must be always t_end-t_begin +1

	ts(int tb, int te) {
		data = NULL;
		if (tb>te) return;

        t_begin = tb;
        t_end = te;

        srand(time(NULL));
        N = te-tb+1;
        data = new double[N];
        for (int i=0; i<N; i++) data[i] = double(rand())/RAND_MAX;
	}

	~ts() {
		if (data != NULL) delete[] data;
	}

	void print() {
		for (int i=0; i<N; i++) {
			cout << t_begin+i << " " << data[i] << endl;
		}
	}

	double operator[](int t) {
		if (t<t_begin) return 0.0;
		if (t>t_end) return 0.0;
		return data[t-t_begin];
	}


	int length() {
        if (t_end-t_begin>0) return t_end-t_begin+1;
        return 0;
    }

	void moving_average(int L) {
		double* new_data = new double[length()+L-1];
		double* temp = data;
		data = new_data;

		for (int i=0; i<length()+L-1; i++) {
			data[i] = 0.0;
			for (int j=0; j<L; j++) {
				if (i-j<0) break;
				data[i]+=temp[i-j];
			}
			data[i]/=L;
		}

		t_end +=L-1;
		N += L-1;

		delete[] temp;
		return;
	}


}; // end of ts

ts& operator +=(ts& num1, ts& num2)
{
    int start = min(num1.t_begin, num2.t_begin);
    int finish = max(num1.t_end, num2.t_end);

    int N = finish-start+1;
    int n1 = num1.length();
    int n2 = num2.length();


    double* data = new double[N];
    double* p=data;
    for (int t=start; t<=finish; t++,p++) {
    	*p = num1[t] + num2[t];
    }

    num1.t_begin=start;
    num1.t_end=finish;
    num1.N=N;

    delete[] num1.data;

    num1.data = data;

    return num1;
}

int main() {
//	ts a(120,129);
	ts b(123, 128);

//	a.print(); cout&lt;&lt;endl;

	b.print(); cout<<endl;

	b.moving_average(3);
	b.print(); cout<<endl;

//	a+=b;
//	a.print();
    //system("PAUSE");
}

Assignment for 24.11.2015: Extend the ts (time-series) object below with an operator such that we can add two time series by matching their indices. Your program should extend the array if neccessary.

Example
// Suppose a_3 = 1, a_4 = 3, a_5 = 2
ts a(3,5);
// Suppose b_5 = 4, b_6 = 3
ts b(5,6);
// Suppose c_10 = 6
ts c(10,10);


a+=b
// a is now
// a_3 = 1, a_4 = 3, a_5 = 6, a_6 = 3

a+=c
// a is now
// a_3 = 1, a_4 = 3, a_5 = 6, a_6 = 3, a_7 = 0, a_8 = 0, a_9 = 0, a_10 = 6

Programs covered at 14.11.2015

Time Series, Dynamic Memory allocation
#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

struct ts {
	int t_begin;
	int t_end;
	double * data;
	int N; // length

	ts(int tb, int te) {

		data = NULL;
		if (tb>te) return;
		t_begin = tb;
		t_end = te;
		srand(time(NULL));

		N = te-tb+1;
		data = new double[N];

		for (int i=0; i<N; i++) data[i] = double(rand())/RAND_MAX;
	}

	~ts() {
		if (data != NULL) delete[] data;
	}

	void print() {

		for (int i=0; i<N; i++) {
			cout << t_begin+i << " " << data[i] << endl;
		}

	}

};



int main() {
	ts a(120,124);
	ts b(123, 147);
	a.print();
	b.print();

}
Pointers
#include <iostream>
using namespace std;

int main() {
	double x = 2.3;

	double* xp = new double[12];

	for (double* p=xp; p<xp+12; p++) *p = 0;

	for (double* p=xp; p<xp+12; p++) cout << *p << endl;

	for (int i=0; i<12; i++) cout << xp[i] << endl;

	delete[] xp;

	return 0;
}

Assignment for 14.11.2015: Reimplement Sorting and Binary search for a _Pair_ object, defined as a structure (see below). The pair keeps the data and the original location. Given an array of Pair's, your implementation should return the index of the number in the original, unsorted, array (a number between 0 and arraySize). If the numberSearched is not in the array, you return -1.

Programs covered at 07.11.2015

Structures, Constructors, Methods and operator overloading
//============================================================================
// Name        : arrays4.cpp
// Author      : A. Taylan Cemgil
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include <iostream>
using namespace std;

struct Pair {
	double x;
	int i;

	Pair (int ii = 0, double xx = 0) {
		i = ii;
		x = xx;
	}

	void print() {
		cout << "(" << i << "," << x << ")" << endl;
	}

};

ostream& operator<<(ostream& os, Pair& p) {
	os << "(" << p.i << "," << p.x << ")";
	return os;
}


Pair& operator+=(Pair& a, Pair& b) {
	a.i += b.i;
	a.x += b.x;

	return a;
}


int main() {
	int j = 3;
	int k = j<<2;
	cout << k << endl;

	Pair c;
	Pair a(0, 3);
	Pair b(6, 12.4);

	cout << a << endl;

	a += b += c; // ==> operator+=(a, b)

	//a.print();

	//b.print();

	cout << a << endl;
}
Computing and operator overloading with Complex numbers
//============================================================================
// Name        : complex_numbers.cpp
// Author      : A. Taylan Cemgil
//============================================================================

#include <iostream>
#include <cmath>

using namespace std;

struct Complex {
	double re;
	double im;
	Complex(double r=0, double i=0) {
		re = r;
		im = i;
	}

	double magnitude() {
		return sqrt(re*re + im*im);
	}

	double phase() {
		return atan(im/re);
	}

};

Complex operator+(Complex& a, Complex& b) {
	return Complex(a.re+b.re, a.im+b.im);
};

Complex operator*(Complex& a, Complex& b) {
	double r = a.re*b.re - a.im*b.im;
	double i = a.re*b.im + a.im*b.re;
	return Complex(r, i);
};


ostream& operator<<(ostream& os, const Complex c) {
	os << "(" << c.re << "+" << c.im << "j)";
	return os;
}


int main() {
	Complex c, d;
	Complex a(2, 3);
	Complex b(-2, 3);

	c = a+b;
	d = a*b;

	cout << (a+b) << endl;
	cout << c << " " << c.magnitude() << " " << c.phase() << endl;
	cout << d << " " << d.magnitude() << " " << d.phase() << endl;
}

Assignment for 7.11.2015: Implement Binary search. Given an array of numbers, your implementation should return the index of the number in the array (a number between 0 and arraySize). If the numberSearched is not in the array, you return -1.

index = binary_search(x, arraySize, numberSearched)

Programs covered at 31.10.2015

Computing various statistics of an array of numbers
//============================================================================
// Name        : arrays3.cpp
// Author      : A. Taylan Cemgil
// Description : Computes various statistics given an array of numbers
//============================================================================

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>

// using namespace std;

unsigned int N = 11;
unsigned int MAX_NUM = 1000;

void print_array(double x[], unsigned int N) {
	for (unsigned int i=0; i<N; i++) {
		std::cout << "(" << i<< "," << x[i] << ")";
	}
	std::cout << std::endl;
}

double max(double x[], unsigned int N) {
	double mx = 0;
	for (int i=0; i<N; i++) {
		if (x[i]>mx) mx = x[i];
	}
	return mx;
}

double min(double x[], unsigned int N) {
	double mn = MAX_NUM;
	for (unsigned int i=0; i<N; i++) {
		if (x[i]<mn) mn = x[i];
	}
	return mn;
}

double mean(double x[], unsigned int N) {
	double mu = 0;
	for (unsigned int i=0; i<N; i++) {
		mu = x[i]/(i+1) + i*mu/(i+1);
	}
	return mu;
}

double stddev(double x[], unsigned int N) {
	double mu = 0;
	double ss = 0;
	for (unsigned int i=0; i<N; i++) {
		mu = x[i]/(i+1) + i*mu/(i+1);
		ss = x[i]*x[i]/(i+1) + i*ss/(i+1);
	}
	return sqrt(ss - mu*mu);
}

// Insertion sort
void sort(double x[], unsigned int N) {
	for (unsigned int i=1;i<N;i++) {
		double d=x[i];
		int j;
		for (j = i-1 ; j>=0; j--) {
			if (d>=x[j]) break;
			x[j+1] = x[j];
		}
		x[j+1] = d;
	}
}


int main() {
	double x[N];

//	srand(time(NULL));
	srand(0);


	for (unsigned int i=0; i<N; i++) {
			x[i] = rand() % MAX_NUM;
	}

	print_array(x, N);
	sort(x, N);
	print_array(x, N);



	std::cout << "Min    : " << min(x, N) << std::endl;
	std::cout << "Max    : " << max(x, N) << std::endl;
	std::cout << "Average: " << mean(x, N) << std::endl;
	std::cout << "Std    : " << stddev(x, N) << std::endl;
	std::cout << "Median    : " << x[N/2] << std::endl;


	return 0;
}

Assignment for 31.10.2015: Write a program for pricing European Call Options using Monte Carlo assign01.pdf

Programs covered at 24.10.2015

Hitting Probability of a Gaussian random walk

#include <iostream>
#include <cstdlib>  // rand()
#include <ctime>    // time()
#include <cmath>    // cos(), sin(), log()

using namespace std;

double randn() {
	double u1 = double(rand())/RAND_MAX;
	double u2 = double(rand())/RAND_MAX;

	double x1 = -2*log(u1);
	double x2 = 2*M_PI*u2;

	double y1 = sqrt(x1)*cos(x2);
	double y2 = sqrt(x1)*sin(x2);
	return y1;
}

double estimate_ruin_probability(
		double S_0,
		double S_T,
		double T = 1.0,
		double steps = 100,
		double s = 1,
		unsigned long N = 10000
		) {
	// ---------------------------------------

	unsigned long num_of_hits = 0;

	double Delta = T/steps;

	for (unsigned long  i=0; i<N; i++) {
		double x = S_0;
		for (int t=0; t<steps; t++) {
			if (x>=S_T) {
				num_of_hits++;
				break;
			}

			double epsilon = randn()*sqrt(Delta*s);
			x = x + epsilon;
		}
	}

	return double(num_of_hits)/N;

}

int main() {
	srand(time(NULL));

	double T = 25;
	double steps = 250;
	double S_0 = 25;
	double S_T = 30;
	double s = 1;

	//<span class="builtin">double</span> p = estimate_ruin_probability(S_0, S_T, T, steps, s);
	double p2 = estimate_ruin_probability(10, 12, 2);

	// cout &lt;&lt; p &lt;&lt; endl;
	cout << p2 << endl;

	return 0;
}

Programs covered at 17.10.2015

Gamblers Ruin Simulation
//============================================================================
// Name        : random_example.cpp
// Author      :
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include <iostream>
#include <cstdlib>  // rand()
#include <ctime>    // time()

using namespace std;

int main() {
	// Set seed
	srand(time(NULL));

	long unsigned int NUMBER_OF_TRIALS = 1000000;
	long unsigned int ruin_count = 0;

	// Time Horizon
	int T = 50;

	// winning probability
	double p = 0.48;

	//cout &lt;&lt; RAND_MAX &lt;&lt; endl;

	double starting_capital = 10;

	for (unsigned long n=0; n<NUMBER_OF_TRIALS; n++) {
		int N;
		double capital = starting_capital;
		for (N=0; N<T; N++) {
			double step = double(rand())/RAND_MAX < p ? 1.0 : -1.0;
			capital = capital + step;
			if (capital <= 0) {
				ruin_count++;
				break;
			}
		}
	}

	cout << "Ruin Probability " << double(ruin_count)/NUMBER_OF_TRIALS << endl;

}
Random Coins Generation
#include <iostream>  // for cout
#include <cstdlib>   // for rand() and srand()
#include <ctime>     // for time, randomize the seed

using namespace std;

int main() {

	srand(time(NULL));

	cout << "RAND_MAX = " << RAND_MAX << endl;
	for (int i=0; i<50; i++) {
		cout << rand()%2;
	}

	return 0;
}

Elementary Examples

Example 1
#include <iostream>

using namespace std;

int square(int x0_input) {
	return x0_input*x0_input;
};

int cube(int x0_input) {
	return x0_input*x0_input*x0_input;
};

int main ()
{
	int i = 0;
	double d = 3.14;
	char ch = 'a';

	cout << "Goodby Moon" << " " << square(5) << " " << cube(5) << endl;
	cout << i << endl;
	cout << d << endl;
	cout << ch << endl;

	return 0;
};
Example 2
#include <iostream>
using namespace std;

int main ()
{
	char c;
	for (int i=1; i<=10; i++) {
		cout << i << " " << i*i << " " << i*i*i << endl;
	};
	cin >> c;
	return 0;
};
Example 3
#include <iostream>
using namespace std;

// Calculate f_{n} = sum_i=1^n i and prints <span class="statement">for</span> i=1..n
int main ()
{
	int n = 20;
	int f;

	f = 0;
	for (int i=1; i<=n; i++) {
		f = f + i;
		cout << i << "," << f << endl;
	};
};
Example 4
#include <iostream>

using namespace std;



// Comments come here

/*

*/

// f_{n} = sum_i=1^n i
// f_1 = 1
// f_2 = 1 + 2
// f_3 = 1 + 2 + 3

// f_1 = 1
// f_2 = f_1 + 2
// f_3 = f_2 + 3

int main ()
{

	for (int i=0; i<=65536; i++) {
		cout << i << endl;
	};


};