Vectorized Monte Carlo for Correlated Basket Option Pricing

# Is this the right way to accelerate my Monte-Carlo Simulation


# Is this the right way to accelerate my Monte-Carlo Simulation












I am trying to develop a pricer for Call VS Call and I'm using MonteCarlo method to do so because my stocks are correlated between each others.

Basically my inputs are `stocks`which is a dataframe whose columns are stock name, volatility, weight for the basket, strike, spot and dividendes, `data`which is a dictionary with the risk free rate, the number of simulations to do and the time to maturity and finally `correlation`which is going to be the implied correlation between stocks (I chose to take it as a float instead of correlation matrix because there is usually only two stocks).

I do Monte-Carlo because there is correlation between my stocks but my option is not path dependent, so I thought I could skip generating X steps for each simulation and directly generate the move over the life of the option. Then I could generate an random array of length (number of stocks, number of simulation) and thus get rid of the usual double loop which is slowing the simulation. When I do like this, I'm able to generate about a million simulation in around 1 second.

However, I struggle to really understand if my way of doing it is good and if this approach is theoretically valid.

Here is my code so each of you can see how I do it :

```
import random 
import matplotlib as plt 
import numpy as np 
from math import * 
import pandas as pd 

class CallVsCallCalculation() : 
    def __init__(self,stocks,data,correlation):
        self.stocks = stocks 
        self.data = data 
        self.correlation = correlation 

        no_obs = int(self.data["correlation"]) 
        r = self.data["rate"]
        T= self.data["time_to_maturity"]
        self.weight = np.array(self.stocks["Weight"]).reshape(len(self.stocks["Weight"]),1)
        self.strike= np.array(self.stocks["Strike"])
        self.spot =np.array(self.stocks["Spot"])
        self.dividend =np.array(self.stocks["Dividendes"])
        means= [0]*len(self.stocks["Variance"])
        no_cols=len(self.stocks["Asset_name"])
        self.vol=list(self.stocks["Variance"])
        sd = np.diag(self.vol)
                  

        observations=np.random.normal(0, 1,(no_cols, no_obs))

        cor_matrix = np.zeros(shape=(no_cols,no_cols))
        for i in range (no_cols):
            for j in range (no_cols):
               if i==j :
                    cor_matrix[i][j]=11
                else :
                    cor_matrix[i][j])=self.correlation

        cov_matrix = np.dot(sd,np.dot(cor_matrix, sd))
        Chol = np.linalg.cholesky(cov_matrix)
        self.chol = np.linalg.cholesky(cov_matrix)
        sam_eq_mean = chol.dot(observations)
        s = sam_eq_mean.transpose() + means
        samples = s.transpose()

        std = np.array(self.vol)
        std **=2
        std = -std/2
        std += r
        std -= self.dividend
        std*=(T)

        eps_array = samples*np.sqrt(T) 
        eps_array = eps_array*np.array(self.vol).reshape((2,1))
        expi = std.reshape(len(self.stocks["Variance"]),1) + eps_array
        #Computing the last last price, or the price at expiration
        self.last_price =  self.spot.reshape(len(self.stocks["Variance"]),1)*np.exp(expi) 
        #Dividing each price by the initial price
        self.last_price_perc = self.last_price/self.spot.reshape(len(self.stocks["Variance"]),1)
        #Getting each price regarding its weights
        self.last_price_weighted = self.weight.T @self.last_price_perc 
        #Calculating the strike
        self.strike_basket = (self.weight.T@self.strike)/((self.weight.T@self.spot)[0]) 
        #Discounting the striketo the basket
        self.basket_minus_strike = self.last_price_weighted-self.strike_basket
        #Getting the payoff for the basket
        self.basket_payoff = np.maximum(self.basket_minus_strike,e)
        self.short_payoff = self.basket_payoff.mean()
                            
        #Getting the strike in term of percentage
        self.strike_stock - self.strike/self.spot
        #weighing the last price by the initial price
        self.last_price_on_spot = self.last_price/self.spot.reshape(len(self.stocks["Variance"]),1) 
        #Discounting the strike
        self.last_price_minus_strike = self.last_price_on_spot-self.strike_stock.reshape(len(self.stocks["Variance"]),1)
        #Getting the payoff
        self.payoff_long = np.maximum(self.last_price_minus_strike,0)
        #Taking the mean of each stocks
        self.payoff_long_mean = self.payoff_long.mean(axis=1)
        #Getting also the weight
        self.long_payoff = self.weight.T@self.payoff_long_mean
                                
        #Getting the call price
        self.call_price = (self.long_payoff-self.short_payoff)[0]*np.exp(-r*T)
```

I also get used of that topic to ask a smaller question about the stock price formula :

I struggle to understand what is the real origin of W in this formula when you consider correlated basket. Is it only the random samples and is it random samples that were multiplied by the lower matrix of the Cholesky decomposition of the covariance matrix ?

Thanks everyone for your help on both question, I hope I was clear enough !