programming - 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 
#include 
#include 
#include 
#include 

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 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.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[0] = 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.

Answer

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.

  
  


• 
  

  Consider your code fragment

  
  
  

  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