Información de fondo:
El algoritmo crudo de Monte Carlo para la opción de compra aritmética asiática es $$Y = e^{-rT}(\overline{S}_A - K)^{+}$$ y el control es $$C e^{-rT}(\overline{S}_G - K)^{+}$$ El estimador de la variable de control es
\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*}
Estimamos que $\mathbb{E}[Y(\beta)]$ utilizando la media de la muestra $$\frac{1}{N}\sum Y_i - \hat{\beta}_N^{*}(C_i - \mathbb{E}[C])$$ donde $$\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}}$$
La opción aritmética de Asain con seguimiento discreto tiene funciones de pago \begin{align*} \text{call payoff} &: \ (\overline{S}_A - K)^{+}\\ \text{put payoff} &: \ (K - \overline{S}_A)^{+} \end{align*} donde $$\overline{S}_A = \frac{1}{n}\sum_{i=1}^{n}S(t_i)$$
Las opciones geométricas de Asain tienen las mismas funciones de retribución pero la media aritmética $\overline{S}_A$ se sustituye por la media geométrica $$\overline{S}_G = \left(\prod_{i=1}^{n}S(t_i)\right)^{1/n}$$ El modelo de precios de las acciones es el modelo GBM: $$S(t_i) = S(0)\exp((r-\sigma^2/2)t_i + \sigma W(t_i))$$ donde $W(t)$ es el movimiento browniano estándar en $[0,T]$ . Multiplicando el modelo GBM como $i = 1,\ldots,n$ obtenemos $$\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)$$
Pregunta:
Consideremos una opción de compra aritmética de Asain con $K = 100$ , $r = 0.2$ , $S_0 = 100$ , $\sigma = 0.4$ utilizando el modelo lognormal para el proceso subyacente de las acciones. Dejemos que el vencimiento $T = 0.1$ con $n = 10$ precios en la media. Generar $40$ estimaciones del precio de la opción, utilizando $N = 100,1000,10000$ trayectorias de los precios, con el crudo Monte Carlo.
Mi pregunta con respecto a este problema es la siguiente. Sé que cuando genero mi vector de tiempo $t$ va a ser de tamaño $10$ . Ahora, mi confusión es con respecto al modelo de precio de las acciones $$S(t_i) = S(0)\exp((r-\sigma^2/2)t_i + \sigma W(t_i))$$ ¿Cuál es el tamaño $S(t_i)$ ? Queremos generar $N = 100,1000,10000$ rutas de precios pero mi vector de tiempo $t$ es sólo de tamaño $10$ entonces, ¿cuál será el tamaño de $S(t_i)$ ? Tengo la misma pregunta para $W(t_i)$ . Cualquier sugerencia sobre lo siguiente será muy apreciada.
Actualización:
Creo que el tamaño de $S(t_i)$ y $W(t_i)$ será de tamaño $10$ pero no entiendo cómo vamos a generar el $N$ caminos de tamaño el tamaño de $N$ es mucho mayor que $10$ .
El modelo lognormal es $$S(t_n) = S(0)\exp\left((r-\sigma^2/2)t_n \sigma\sqrt{t_n}Z\right)$$
No estoy seguro de cómo todavía para generar esto por lo que el tamaño que necesito para establecer $Z$ ...etc...
Este es el código que he hecho hasta ahora, sería genial si alguien puede decirme cómo genero el $40$ estimaciones del precio de la opción utilizando $N = 100,1000,10000$ rutas de precios.
// 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 pricesEste es mi código completo, he podido resolverlo yo mismo:
#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();
}
}
