# Use random-shift Halton sequence to obtain 40 independent estimates for the price of a European call

Background Information:

Random-shift Halton sequence: Consider the first six Halton vectors in dimension $2$, using base $2$ and $3$:

$$\begin{bmatrix} 1/2\\ 1/3 \end{bmatrix}, \begin{bmatrix} 1/4\\ 2/3 \end{bmatrix},\begin{bmatrix} 3/4\\ 1/9 \end{bmatrix},\begin{bmatrix} 1/8\\ 4/9 \end{bmatrix},\begin{bmatrix} 5/8\\ 7/9 \end{bmatrix}, \begin{bmatrix} 3/8\\ 2/9 \end{bmatrix}$$

To obtain independent randomization of this set, generate two random numbers $u_1$,$u_2$ from $U(0,1)$, say, $u_1 = 0.7$ and $u_2 = 0.1$. Add the vector $u = \begin{bmatrix} u_1\\ u_2 \end{bmatrix}$ to each vector above, componentwise, and take the fractional part of the sum. For example, the first vector becomes $$\begin{bmatrix} \{1/2 + 0.7\}\\ \{1/3 + 0.1\} \end{bmatrix}= \begin{bmatrix} 0.2\\ 0.4\overline{3} \end{bmatrix}$$ (here $\{x\}$ denotes the fractional part of $x$). Similarly, add $u$ to the rest of the vectors. Here is the Pseudocode for simulating the geometric Brownian motion model: Question:

Consider a European call option with parameters: $T = 1$, $K = 50$, $r = 0.1$, $\sigma = 0.3$, and $S_0 = 50$

Use random-shift Halton sequences to obtain $40$ "independent" estimates for the price of the option. For each estimate, use $N = 10,000$ price paths. To simulate a path, you will simulate the geometric Brownian motion model with $\mu = r$, and using $10$ time steps $t_0 = 0$, $t_1 = \Delta t$, $t_2 = 2\Delta t,\ldots,t_{10} = 10\Delta t = T$. Use the Box-Muller method to generate the standard normal numbers.

I have been able to sucessfully use the Box-Muller method to generate the standard normal numbers which I checked are correct. I also have simulate the geometric Brownian motion model with $\mu = r$ but I am not completely sure it is correct.

The part that I am confused about is assuming the geometric Brownian motion model is correct how I apply the random-shift Halton sequences to obtain the $40$ "independent" estimates for the price of the option.

I programmed the latter in C++. Here is the code which should be pretty easy to read:

#include <iostream>
#include <cmath>
#include <math.h>
#include <fstream>
#include <random>

using namespace std;

double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(int n, int N,double u_1, double u_2);
double Box_MullerY(int n, int N,double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);

int main(){
int n = 10;
int N = 10000;

// Calculate BS
double T = 1.0;
double K = 50.0;
double r = 0.1;
double sigma = 0.3;
double S0 = 50.0;
double price = Black_Scholes(T,K,S0,r,sigma);
//cout << price << endl;

// Random-shift Halton Sequence
Random_Shift_Halton(n,N,r,sigma,T, S0);

}

double phi(double x){
return 0.5 * erfc(-x * M_SQRT1_2);
}

double Black_Scholes(double T, double K, double S0, double r, double sigma){
double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
return C;
}

double Halton_Seq(int index, double base){
double f = 1.0, r = 0.0;
while(index > 0){
f = f/base;
r = r + f*(fmod(index,base));
index = index/base;
}
return r;
}

double Box_MullerX(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double X = sqrt(R)*cos(theta);

return X;
}

double Box_MullerY(int n, int N,double u_1, double u_2){
double R = -2.0*log(u_1);
double theta = 2*M_PI*u_2;
double Y = sqrt(R)*sin(theta);

return Y;
}

double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
//  Generate the Halton_Sequence
double hs[N + 1];
for(int i = 0; i <= N; i++){
hs[i] = Halton_Seq(i,2.0);
}
// Generate two random numbers from U(0,1)
double u[N+1];
random_device rd;
mt19937 e2(rd());
uniform_real_distribution<double> dist(0, 1);
for(int i = 0; i <= N; i++){
u[i] = dist(e2);
}
// Add the vector u to each vector hs and take fraction part of the sum
double shifted_hs[N+1];
for(int i = 0; i <= N; i++){
shifted_hs[i] = hs[i] + u[i];
if(shifted_hs[i] > 1){
shifted_hs[i] = shifted_hs[i] - 1;
}
}

// Use Geometric Brownian Motion
double Z[N];
for(int i = 0; i < N; i++){
if(i % 2 == 0){
Z[i] = Box_MullerX(n,N,shifted_hs[i+1],shifted_hs[i+2]);
}else{
Z[i] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
}
}
double ZZ[N + 1];
ZZ = 0.0;
for(int i = 1; i <= N; i++){
ZZ[i] = Z[i-1];
}

for(int i = 1; i <= N; i++){
//cout << ZZ[i] << endl;
}

double S[N + 1];
S = S0;
double t[n+1];
for(int i = 0; i <= n; i++){
t[i] = (double)i*T/(n-1);
}
for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[i] = S[i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*ZZ[i]);
}
}

for(int i = 1; i <= n; i++){
//cout << S[i] << endl;
}

// Use random-shift halton sequence to obtain 40 independent estimates for the price of the option

}

}


If any one has any suggestions on what I do from here I would greatly appreciate it.

• Is this homework/an assignment? I am asking because apart from your question, there are multiple things wrong with your other code. – LocalVolatility Mar 11 '17 at 4:52

Actual Question

I suppose this is homework, so I will only outline the steps. The way I understand this question is as follows:

1. Build a function that simulates different 10,000 sample paths of your underlying asset with 10 equidistant time steps each.
2. For each of these paths, compute the terminal European call option payoff. Their discounted sample mean is your estimate for the European call price.

Repeat steps 1 and 2 for a total of 40 times with different randomized offset for the Halton sequences to get 40 different estimates and record the results.

Other Code

Looking at your code, I can see quite a few issues. Here are some hints to get you started:

• You don't carry out the randomization of the Halton sequence at all. This should happen between generating standard Halton sequences and using them in the Box-Muller transform to generate random normals. The excerpt you posted is very clear on how it is done.

• This code fragment

for(int i = 0; i < N; i++){
if(i % 2 == 0){
Z[i] = Box_MullerX(n,N,hs[i+1],hs[i+2]);
}else{
Z[i] = Box_MullerY(n,N,hs[i],hs[i+1]);
}
}


only happens to not go out-of-bounds because N is even in your main.

• If you want to generate N sample paths, then you'd want S to have dimension N if you just want to keep track of the current spot of each path (which is enough). I.e. S[j] is the current spot of the $j$-th path and you'd have to initialize all entries in S to S0.

for(int j = 1; j <= N; j++){
for(int i = 1; i <= n; i++){
S[i] = S[i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*ZZ[i]);
}
}


The outer loop does nothing except for heating up your CPU. You never use the loop variable j but simply run the exact same inner loop N times. Put differently, you simulate a single path only instead of N.

• The Box-Muller transform needs independent uniform on the unit square. You use two successive items from the same Halton sequence. If you create a scatterplot of the supposedly normal variates, you get something like the first of the below plots. The second was produced using two Halton sequences with different bases. See also the corresponding discussion in Chapter 9.3.1 of Jaeckel (2002).

• References

Jaeckel, Peter (2002) Monte Carlo Methods in Financial Engineering: John Wiley & Sons

• Thank you for outlining it for me. Yes, it is homework. I have some questions on your comments. You stated that I don't do the randomization of the Halton sequence. Are you saying that when I generate the Halton sequence I need to store them in a randomized way? Could you clarify what you mean about tha – Wolfy Mar 11 '17 at 6:10
• ** Could you clarify what you mean by the first code block you mentioned. I was just trying to store he $X$'s and $Y$'s in order into the array $Z$. – Wolfy Mar 11 '17 at 6:17
• 1) In your own description it says "add the vector $u$ to each vector above" - you simply do not do this. This has to happen before the Box-Muller transform. Note that this is yet a different issue than the one in my last point. 2) Check what indices of hs you are accessing when N is odd - let's say N = 3. You'll find that hs has size 4 and you try to access the element with index 4 which is out of range. – LocalVolatility Mar 11 '17 at 11:12
• Ok, well $N = 10,000$ so it should be fine for that case since $N$ will not be odd ever. – Wolfy Mar 11 '17 at 16:28
• So to clarify, I have to add the vector $u$ to each vector above before I apply the Box-Mueller tranform, then I proceed to simulating the geometric Brownian motion? – Wolfy Mar 11 '17 at 16:32