Crude Monte Carlo Simulation for Arithmetic Asian Options

# Using crude Monte Carlo


# Using crude Monte Carlo












Background Information:

The crude Monte Carlo algorithm for the arithmetic Asian call option is $$Y = e^{-rT}(\overline{S}_A - K)^{+}$$ and the control is $$C e^{-rT}(\overline{S}_G - K)^{+}$$ The control variate estimator is

\begin{align*} Y(\beta) &= Y - \beta(C - \mu_C)\\ &= e^{-rT}(\overline{S}_A - K)^{+} - \beta(e^{-rT}(\overline{S}_G - K)^{+} - \mathbb{E}[C]) \end{align*}

We estimate $\mathbb{E}[Y(\beta)]$ using the sample mean $$\frac{1}{N}\sum Y_i - \hat{\beta}_N^{*}(C_i - \mathbb{E}[C])$$ where $$\hat{\beta}_N^{*} = \frac{\sum_{i=1}^{n}(Y_i - \overline{Y})(C_i - \overline{C})}{\sum_{i=1}^{n}(C_i - \overline{C})^{2}}$$

Arithmetic Asain option with discrete monitoring have payoff functions \begin{align*} \text{call payoff} &: \ (\overline{S}_A - K)^{+}\\ \text{put payoff} &: \ (K - \overline{S}_A)^{+} \end{align*} where $$\overline{S}_A = \frac{1}{n}\sum_{i=1}^{n}S(t_i)$$

Geometric Asain options have the same payoff functions but the arithmetic Asain options have the same payoff function but the arithmetic average $\overline{S}_A$ is replaced by the geometric average $$\overline{S}_G = \left(\prod_{i=1}^{n}S(t_i)\right)^{1/n}$$ The stock price model is the GBM model: $$S(t_i) = S(0)\exp((r-\sigma^2/2)t_i + \sigma W(t_i))$$ where $W(t)$ is the standard Brownian motion on $[0,T]$. Multiplying the GBM model as $i = 1,\ldots,n$ we obtain $$\overline{S}_G = \left(\prod_{i=1}^{n}S(t_i)\right)^{1/n} = S(0)\exp\left(\left(r - \frac{\sigma^2}{2}\right)\left(\frac{1}{n}\sum_{i=1}^{n}t_i\right) + \frac{\sigma}{n}\sum_{i=1}^{n}W(t_i)\right)$$

Question:

> Consider an Asain arithmetic call option with $K = 100$, $r = 0.2$, $S_0 = 100$, $\sigma = 0.4$, using the lognormal model for the underlying stock process. Let the expiry $T = 0.1$ with $n = 10$ prices in the average. Generate $40$ estimates for the option price, using $N = 100,1000,10000$ price paths, with crude Monte Carlo.

My question in regards to this problem is the following. I know that when I generate my time vector $t$ it is going to be of size $10$. Now, my confusion is in regards to the stock price model $$S(t_i) = S(0)\exp((r-\sigma^2/2)t_i + \sigma W(t_i))$$ What is the size $S(t_i)$? We want to generate $N = 100,1000,10000$ price paths but my time vector $t$ is only of size $10$ so what will be the size of $S(t_i)$? I have the same question for $W(t_i)$. Any suggestions on the following are greatly appreciated.

Update:

I believe that the size of $S(t_i)$ and $W(t_i)$ will be of size $10$ but I do not understand how we will generate the $N$ paths size the size of $N$ is much larger than $10$.

The lognormal model is $$S(t_n) = S(0)\exp\left((r-\sigma^2/2)t_n \sigma\sqrt{t_n}Z\right)$$

I am not sure how to still to generate this because of what size I need to set $Z$, etc...

Here is the code I have done so far, it would be great if someone can just tell me how I generate the $40$ estimates for the option price using $N = 100,1000,10000$ price paths.

```
// Create time vector
VectorXd tt = time_vector(0.0,T,n);
VectorXd t(n);
for(int i = 0; i < n; i++){
        t(i) = tt(i+1);
}

// Generate standard normal Z matrix
MatrixXd Z = generateGaussianNoise(N,n);

// Generate N paths of stock prices
```

This is my complete code, I was able to figure it out myself:

```
#include <iostream>
#include <cmath>
#include <math.h>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <random>
#include <time.h>

using namespace Eigen;
using namespace std;

void crudeMonteCarlo(int N,double K, double r, double S0, double sigma, double T, int n);
VectorXd time_vector(double min, double max, int n);

int main(){
    int N = 100;
    double K = 100;
    double r = 0.2;
    double S0 = 100;
    double sigma = 0.4;
    double T = 0.1;
    int n = 10;

    crudeMonteCarlo(N,K,r,S0,sigma,T,n);

    return 0;
}

VectorXd time_vector(double min, double max, int n){
    VectorXd m(n + 1);
     double delta = (max-min)/n;
     for(int i = 0; i <= n; i++){
             m(i) = min + i*delta;
     }
    return m;
}

MatrixXd generateGaussianNoise(int M, int N){
    MatrixXd Z(M,N);
    static random_device rd;
    static mt19937 e2(time(0));
    normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        for(int j = 0; j < N; j++){
            Z(i,j) = dist(e2);
        }
    }
    return Z;
}

void crudeMonteCarlo(int N,double K, double r, double S0, double sigma, double T, int n){
    // Create time vector
    VectorXd tt = time_vector(0.0,T,n);
    VectorXd t(n);
    double dt = T/n;
    for(int i = 0; i < n; i++){
            t(i) = tt(i+1);
    }

    // Generate standard normal Z matrix
    //MatrixXd Z = generateGaussianNoise(N,n);

    // Generate 40 estimates for the option price, using N paths
    int m = 40;
    MatrixXd SS(N,n+1);
    VectorXd S(m);
    for(int k = 0; k < m; k++){
        MatrixXd Z = generateGaussianNoise(N,n);
        for(int i = 0; i < N; i++){
            SS(i,0) = S0;
            for(int j = 1; j <= n; j++){
                SS(i,j) = SS(i,j-1)*exp((double) (r - pow(sigma,2.0))*dt + sigma*sqrt(dt)*(double)Z(i,j-1));
            }
        }
        S(k) = SS.mean();
    }

}
```

## Answer by will (score 4, accepted)

https://quant.stackexchange.com/a/33461

Here's some pseudo code to generate your valuations:

```
Asset asset(); //some asset object that holds all the data.

vector<double> timeline; //this is your vector of times.
size_t nTimes = timeline.size();
vector<double> dts(nTimes); //vector of fwd prices on the same timeline.
vector<double> fwds(nTimes); //vector of fwd prices on the same timeline.
vector<double> fwdRatios(nTimes); //vector of fwd prices on the same timeline.
double vol = 0.2; //Your BS vol.    

for(size_t iTimeline=0; iTimeline<nTimes; iTimeline++){
    fwds(i) = asset.getForward(timeline[iTimeline]);
    if(iTimeline==0){
        fwdRatios(0) = fwds(i) / asset.getSpot();
        dts(0) = timeline(0);
    }else{
        fwdRatios(i) = fwds(i) / fwds(i-1);
        dts(i) = timeline(i) - timeline(i-1);
    }
}

//Some random number generator - 
PRNG<Distribution::Normal> prng();

size_t nTrajectories = 100000;

Matrix<double> trajectories(nTimes, nTrajectories);
for(size_t iTrajectory=0; iTrajectory<nTrajectories ; iTrajectory++){
    trajectories(0, iTrajectory) = asset.getSpot() * fwdRatios(0) * exp(-0.5*vol*vol*dts(0) + vol*sqrt(dts(0)*prng.next());

    for(size_t iTimeline=1; iTimeline<nTimes; iTimeline++){
        trajectories(iTimeline, iTrajectory) = trajectories(iTimeline-1, iTrajectory) * fwdRatios(iTimeline) * exp(-0.5*vol*vol*dts(iTimeline) + vol*sqrt(dts(iTimeline))*prng.next());
    }
}

vector<double> payoffs(nTrajectories);
for(size_t iTrajectory=0; iTrajectory<nTrajectories ; iTrajectory++){
    payoffs(iTrajectory) = calculatePayoff(trajectories.row(iTrajectory));
}
```

That'll give you one sample price. If you want to do it 40 times, just put another foor loop around it to do it 40 times.

You'll not be able to copy paste that and compile it though, because i just typed it out and made up functions...