Sampling from Cauchy (and Gaussian) distribution

How we construct  a random draw  from a Cauchy distribution, if we start with a draw from uniform distribution? That’s what asks the Exercise 7.3 of the book The Concepts and Practice of Mathematical Finance by Mark Joshi. We post here some C++ code as an answer. This is easily done by explicitly computing the inverse cumulative distribution function F^{-1}(x) of the Cauchy distribution and plugging then random numbers generated with a uniform distribution (this is done with the function GetOneRandomByICD).

f(x) = \frac{1}{\pi}\frac{1}{1+x^2},\; F(x) = \frac{1}{\pi}\arctan (x) + \frac{1}{2},
F^{-1}(x) = \tan [\pi (x-1/2)]

More tricky is the case in which we don’t have the F^{-1}(x). In that case we can obtain our sample by using two random uniforms – GetOneCauchyByUniforms –  or using two random numbers with gaussian distribution – GetOneCauchyByGaussians. In the latter case of course we need also a routine  to compute gaussian numbers. We present here the routine GetOneGaussianByRatio, which is similar to GetOneCauchyByUniforms (both are based on ratio of uniforms – see Probability and Random Processes by Grimmet and Stirzaker). The mathematical proofs are in the probability section of my notes.

//Random.h
#ifndef RANDOM_H
#define RANDOM_H

double CauchyICD(double uniform);
double GetOneRandomByICD(double (*func_ptr)(double));

double GetOneCauchyByUniforms();
double GetOneCauchyByGaussian();

double GetOneGaussianByRatio();

#endif
//Random.cpp
#include "Random.h"
#include <cmath>
#include <cstdlib>

double CauchyICD(double x){
 return std::tan(M_PI*(x-0.5));
}
double GetOneRandomByICD(double (*func_ptr)(double)){
	return (*func_ptr)(rand()/static_cast<double>(RAND_MAX));
	}	

double GetOneGaussianByRatio(){

	double x;
	double y;

	double ratio;

	do {
		x = 2.*rand()/static_cast(RAND_MAX)-1.;
		y = 2.*rand()/static_cast(RAND_MAX)-1.;
		ratio = x/y;
		}
	while((y > std::exp(-ratio*ratio*0.25))     
           || (y < -std::exp(-ratio*ratio*0.25)) );

	return ratio;
	}

double GetOneCauchyByUniforms(){

	double x;
	double y;

	double sizeSquared;

	do {
		x = 2.*rand()/static_cast<double>(RAND_MAX)-1.;
		y = 2.*rand()/static_cast<double>(RAND_MAX)-1.;
		sizeSquared = x*x + y*y;
		}
	while(sizeSquared >= 4.0);

	return y/x;
	}

double GetOneCauchyByGaussian(){
	double x = GetOneGaussianByRatio();
	double y = GetOneGaussianByRatio();

	return x/y;
	}

A simple test program which write to files samples of Cauchy distributed numbers:

//testRandom.cpp
#include "Random.h"
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>

int main(){

	long sampleLength = 10000;

	std::ofstream myfile1;
	myfile1.open ("inverse_cauchy.txt");
	std::ofstream myfile2;
	myfile2.open ("ratio_cauchy.txt");	
        std::ofstream myfile3;
	myfile3.open ("gauss_cauchy.txt");	

    for(long n=0; n< sampleLength; n++){
		 myfile1 << GetOneRandomByICD(&CauchyICD) << "\n";
		 myfile2 << GetOneCauchyByRatio() << "\n";
		 myfile3 << GetOneCauchyByGaussian() << "\n";
		}

	myfile1.close();
	myfile2.close();
	myfile3.close();

	return 0;
	}

Now just compile and try it.

Advertisements

Pricing on Trees notes

Just added to the Notes page a few pages about option pricing on trees. We explicitly show the equivalence of replicating, hedging or using a risk-neutral approach to price options on a binomial tree. A Mathematica implementation of binomial and trinomial tree will be added soon, with discussion on convergence issues.

European Options on Google Play

Just uploaded to the market European Options:

Easily compute price and greeks of european-style option derivatives. European Options is an handy app to compute price and greeks of Vanilla, Digital and Barrier european-style options (for definitions we refer to the standard book “Hull – Options, futures and other derivatives”). It is planned to add more exotic options. Price and greeks are computed via analytical formulas.

The source code is available at my GitHub repo.

Get it on Google Play

European options App – source code

Pushed to my GitHub repo a first version of European Options: an handy Android app which computes price and greeks for Vanilla, Digital and Barrier european-style options. Basically it is the Java translation of some of the C++ code I wrote this summer following Joshi books. In a few days I will upload to the market the app.

Hello world!

I’m a PhD in theoretical physics (as one can easily guess from the title of the blog I worked on Higgs physics) and I would like to move towards quantitative finance.

I’ll post here notes, code snippets,  tutorials, book reviews and web resources which I’m collecting while studying quant finance.  In the hope that the resources within this web page will help people with similar background.