*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 of the Cauchy distribution and plugging then random numbers generated with a uniform distribution (this is done with the function

*GetOneRandomByICD*).

More tricky is the case in which we don’t have the . 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.