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