# 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();
}

}

• $S(t_i)$ is of sizee $N$ - but where you have $t_i, i \in 0 \ldots n$, then you have n vectors of price. Typically you'd describe them as paths though, so you have to generate a $N$ paths of prices of length $n$. – will Apr 1 '17 at 19:09
• @will I am still a bit lost by that, could you answer this question providing pseduocode of what you mean? I am writing this in C++. – Wolfy Apr 1 '17 at 19:11
• There are hundreds, if not thousands of black sholes montecarlo implementations on the internet already... for each path, you essentially want to do this: $S_i = S_{i-1} \frac{F_i}{F_{i-1}} e^{-0.5 \sigma^2 \mathrm{d}t + \sigma \sqrt{\mathrm{d}t} X}$ where $X$ is a normally distributed random number. – will Apr 1 '17 at 19:17
• @will and does $i = 1,\ldots,N$ in my question. Also, I am trying to follow what the notes do so if you see the $W(t_i)$ in the question is the standard Brownian motion. Also, instead of $dt$ I have my time vector $t$, also in the question. This is why I am confused, I need to generate $N$ paths but by the equation presented (GBM) I have to use my $t_i$ value which is of size $10$ but I need to create $N$ paths, do you understand how I am confused? – Wolfy Apr 1 '17 at 19:19
• Though it doesn't add much to @will 's comments: The first part of my answer to this (your) question quant.stackexchange.com/questions/32942 pretty much answers your current question as well. By the way - your actual question has nothing to do with Asian options. By adding all those unnecessary details, you make it less clear and will probably attract less answers. As far as I can see, you problem is a very basic one - you seem to confuse "number of steps", "number of paths" and "number of estimates" in Monte Carlo simulations. – LocalVolatility Apr 3 '17 at 17:00

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...