Note: The first version of this file was written in Mathematica 4.2. However, the author hopes that all algorithms and functions were successfully (and regurarly) rewritten to be in line with up-to-date version of Mathematica.

Option pricing via Monte Carlo simulation
(Odhad ceny opce dle simulace Monte Carlo)

Symbols and Notations

Bellow, we review the notations and symbols which are most frequently used in this text.

	spot	spot price of the underlying asset at time zero	spotová cena podkladového aktiva v čase nula
	spotT	spot price of the underlying asset at maturity time	spotová cena podkladového aktiva v době zralosti
K	exe	exercise price (strike price)	realizační cena
T	mat	maturity time	doba zralosti
τ	timeto	time to maturity	doba do zralosti
r	r	riskless rate	bezrizikový výnos
σ	vol	volatility of returns of S	volatilita výnosů S
q	q	dividend yield	dividendový výnos
N	number	number of independent scenarios	počet nezávislých scénářů
□	length	length of the time interval used for simulation	délka časového úseku pro simulaci
θ	driftVG(NIG)	skewness parameter	parametr řídící šikmost
ν	gammavariance	kurtosis parameter	parametr řídící špičatost
ω	omega	mean correcting parameter	parametr upravující střední hodnotu
dist[θ]	distribution	a general notation to a target distribution with parameters θ	obecné označení cílového rozdělení pravděpodobnosti s parametry θ
G[a,b]	gamma dist.	gamma distribution with shape and scale parameters a and b	gama pravděpodobnostní rozdělení s parametry a a b

Required functions of Mathematica: Table, Mean, StandardDeviation, Variance, Skewness, Kurtosis, Module, SeedRandom, RandomReal, Clear, Take, Drop, Append, Delete, Flatten, Timing, Join, Quantile, CharacteristicFunction, PDF, UniformDistribution, Grid, GraphicsArray, Histogram, Plot, NormalDistribution, CDF, GammaDistribution, ...

Some basic facts about probability, random numbers and statistical estimation

This section is under construction...

Random numbers from prespecified distributions

In this section we provide basic moduls to generate random numbers from selected distributions, useful in Financial modeling. However, we start with general moduls, without specifying the target distribution.

Define first the vector of parameters for further utilization (head of resulting tables, head), the function for grid to provide the table of particular results, ResultsGrid, and the parameter to set the maximal number of trials (kk):

General modules

PMC (Plain Monte Carlo)

This is the simplest approach, since we use solely the predefined algorithm of Mathematica specified within the function RandomReal[dist[θ],{number}]. Here, dist[θ] refers to the target distribution with a given set of parameters θ, number is the number of independent scenarios we generate. Within the algorithm bellow, when we specify the internal module for pN, we use nn to allow us to recover exactly the same "random" realizations later, if needed. Thus, the algorithm PMCrannumPar has three variables, dist[θ], nn, and number. It produces the vector of random variables, rannum, for which the estimations of the sample mean, standard deviation, variance, skewness, and kurtosis is obtained. Moreover, the algorithm returns also the number of generated scenarios. Finally, we show how to run the algoritm to get the vector of estimated parameters accompanied by the overall time costs for   random scenarios. Note: We start the definition by clearing the function (PMCrannumPar)

In[39]:=

Alternatively, we can define a simple algorithm that will not produce the estimation of the paramaters, but solely the random numbers. Although it takes more memory to rember the results, it is the only way how to picture (approximate) the probability density function. Basically, we have three choices - the approximation of the PDF by a histogram; smoothing it by a Kernel function; comparing the quantiles (QQ plot). As above, we first realize random scenarios.

In[42]:=

Another way we can proceed is to produce everything within a single algoritm. Subsequently, we pick up separately the parameters of a given random number distribution, including the time costs, and the set of random numbers, see bellow.

To conclude, when we are interested mainly in desriptive statistic parameters, we use the modul with "Par", PMCrannumPar. By contrast, when we need to get the numbers for further utilization, we can use the modul with "Num", PMCrannumNum. Finally, to get both, we use just PMCrannum.

The approach described above can be used in connection with any predefiened distribution. Obviously, if there is a need for a distribution which is not predefined in Mathematica  (eg. VarianceGamma), it must be defined before running the procedure.

AVMC (Monte Carlo with antithetic variates)

The Monte Carlo approach with antithetic variates in its basic form gives an additional path for each original one obtained due to the predefined algorithm of Mathematica, see PMC. Thus, we call only for a given fraction of the required amount of variables, rannum. More particularly, in its simplest form we obtain directly one half of the variables we actually need, rannumpart=number/2 The drawback of the approach is that it can be used directly only for symmetric distributions. If the mean is nonzero, we have to deduce the random realizations from its value MEAN.

Similarly to PMC, we can define two alternative algoritms. First, we state the one to produce the set of random numbers:

And here is the other one:

If the target distribution is not symmetric, we can reformulate the algorithm to produce first random numbers due to the standard uniform distribution, ie we obtain randomly the probabilities, and then apply the inverse transform method.

ITMC (Inverse Transforme Monte Carlo)

Assume that we are able to generate random variables from the standard uniform distribution, but we are interested in some more complex one. In other words, the standard uniform variates provide us randomly with probabilities, ie values of the target distribution function. Assuming that the inverse function to the target distribution function exists, we can transform such variables to obtain the target ones. In Mathematica, we can use the function called Quantile. We state only the basic algorithm.

A strong drawback of ITMC is that the calculation of inverse functions often takes too much time - it is several times more time consuming than PMC.

SSMC (Stratified Sampling Monte Carlo)

The key step of this improvement to simple Monte Carlo approach is to stratify the interval of admissible resutls into subintervals with equal probabilities of occurence (alternativelly, we call these subintervals stratas). In general, there exist two approaches to stratified sampling. We can stratify either the original set or the unit interval, ie probabilites. The latter case is generally more tractable, since it is not required to deal with infinite bounds and the subintervals have equal length. We first define an internal module to provide random numbers due to the standard uniform distribution, pU. Next, we table whole numbers (vector P1) up to the total quantity of scenarios considered, number. Subsequently, we stratify the unit interval and generate one random number for each of them. Finally, we produce random numbers from the target distribution applying function Quantile.

As we will seevlater this algorithm extends ITMC approach to provide more efficient results at almost same time costs. Clearly, the algorithm can be redefined to produce more independent random numbers for each of the stratas. As above, we state the two alternative algorithms bellow. This is the first one:

And here is the other one:

It is obvious, that such approach works only for one-dimensional problems. If more dimensional random vector must be considered, this approach can be easily extended to LHS (Latin Hypercube Sampling) Monte Carlo method. The application is easy. We only need to permute randomly the vector of boundaries for particular stratas P.

AV-SSMC (Stratified Sampling with Antithetic Variates Monte Carlo)

The stratified Monte Carlo approach presented above is based on generating random numbers from Uniform Distribution first. Clearly, we can also think about improving the generator by e.g. AVMC approach. In this subsection we present several alternatives. For simplicity, we will deal only with algorithms of SSMCrannumNum class. We start with algorithm, that provides two paths for each strata: we accompine the original one by a mirrored one, due to the mean of the strata.

Alternatively, both realizations can be symmetric with respect to the mean of the entire interval.

Next choice can be to...

Uniform Distribution

Definition

Uniform distribution indicates that all admissible values can arise with the same probability, ie the probabilities of particular outcomes are equal. In Mathematica, we can get uniformly distributed random numbers calling UniformDistribution[{a, b}], where a and b indicates the bounds of admissible values, b>a. It can be defined either by a characteristic function:

or a distribution function:

Clearly, it is a horizontal line for any x∈[a,b].

Generally we apply a standard uniform distribution, ie uniform distribution with unit margins,  a=0 and  b=1, which we can call for by typing UniformDistribution[{0,1}] or in order to get random variables RandomReal[{}]. Now, we ask Mathematica to calculate analytically the first four moments of the distribution function and we accompany them by the number of scenarios and the time costs.

Algorithms

This part is under construction...

PMC (Plain Monte Carlo)

Now, we utilize the algorithm PMCrannumPar defined above. Before getting the table of results, we must specify kk, the first row head and the last row lr with analytically obtained parameters.

In order to compare the density/distribution we apply the other algorithm.

Now we map the Histogram function to vectors of random numbers generated for particular k.

AVMC (Monte Carlo with antithetic variates)

AVMC approach needs to run only one half of scenarios comparing with PMC. We examine if a decrease in time costs is adequate. We should be also interested, if the estimated parameters converge to the true values rapidly.

Alternatively, we use the other algorithm to compare the overall distribution.

Now we map the Histogram function to vectors of random numbers generated for particular k.

ITMC (Monte Carlo on the basis of inverse transform method)

In order to separate the effect of stratification and the inverse function to the distribution function evaluation, we should run ITMC alone.

SSMC (Stratified Monte Carlo)

Now we can proceed to stratified Monte Carlo. Since it generally shows to be very time consuming to evaluate the inverse function to any distribution function, we switch to the following algorithm:

Now we map the Histogram function to particular vectors of random numbers generated for particular k.

Normal (Gaussian) Distribution

Definition

The term normal distribution is commonly used instead of Gaussian distribution, since such a distribution of random variables is supposed to be the most frequent in practice considerations. In Mathematica, we call for it by typing NormalDistribution[μ, σ], ie with mean μ and its standard deviation σ>0. It can be defined either by a characteristic function:

or a distribution function:

We plot the probability distribution functions for various σ and μ, respectivelly, bellow.

Note: with plain settings the first function is drawn with blue, the second with violet, etc.

Generally we apply a standard normal distribution, ie normal distribution with zero mean and unit variance, for which we can call simply by NormalDistribution[ , ]. The reason is that it can be easilly transformed, with insignificant time costs, into the target distribution by multiplying the original vector by σ and adding μ. Now, we analytically calculate the first four moments of the distribution function and we accompany them by the number of scenarios and the time costs.

Algorithms

This part is under construction...

PMC (Plain Monte Carlo)

Here, we utilize the algorithm PMCrannumPar defined above. Before getting the table of results, we must specify kk, the first row head and the last row lr with analytically obtained parameters.

Next, in order to compare the density/distribution we apply the other algorithm.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

AVMC (Monte Carlo with antithetic variates)

AVMC approach needs to run one half of scenarios comparing with PMC. We examine if a decrease in time costs is adequate. We should be also interested, if the estimated parameters converge to true values rapidly.

Alternatively, we use the other algorithm to compare the overall distribution.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

ITMC (Monte Carlo on the basis of inverse transform method)

In order to separate the effect of stratification and applying the inverse function to the distribution function within SSMC approach, we need to run ITMC alone.

SSMC (Stratified Monte Carlo)

We can proceed now to stratified Monte Carlo. Since it generally seems to be very time consuming to evaluate the inverse function to any distribution function, we switch to the following algorithm:

Now we map the Histogram function to particular vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

Gamma Distribution

Definition

Gamma distribution belongs to the large family of distributions with non-negative realizations. It is therefore useful when modeling of stochastic processes with non-negative realizations (eg. time increments) is considered. In Mathematica, we call for it by GammaDistribution[a, b], ie with parameters of shape a > 0 and scale b > 0. With basic conditions, we can defined it either by a characteristic function:

or a distribution function:

The mean and the variance can be precalculated as follows:

Due to

it follows that and Later, whitin modeling the asset price evolutions, we will require: μ=1 and ν=0.2264. Note also that the Gamma distribution is often defined by parameters α and β with relations to a and b. Under suitable choice of parameters, the Gamma distribution gives us either an Exponential distribution or a Chi-2 distribution.

We plot the probability distribution function bellow.

Note: with plain settings the first function is drawn with blue, the second with violet.

It is important to know that it is sufficient to generate GammaDistribution[a, 1]. The reason is that it can be easilly transformed into target distribution GammaDistribution[a, b] when the results are multiplied by b. Now we analytically calculate the first four moments of the distribution function and we accompany them by the number of scenarios and the time costs, respectivelly.

Algorithms

This part is under construction...

PMC (Plain Monte Carlo)

Here, we utilize the algorithm PMCrannumPar defined above setting the distribution to be GammaDistribution[a,b]. Moreover, before getting the table of results, we must specify kk, the first row head and the last row lr with analytically obtained parameters.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

AVMC (Monte Carlo with antithetic variates)

Since the Gamma probability distribution is not a symmetric one, there is obviously no reason to apply the AVMC method in its basic form.

SSMC (Stratified Monte Carlo)

We can proceed now to stratified Monte Carlo. Since it generally seems to be very time consuming to evaluate the inverse function to any distribution function, we switch directly to the following algorithm to produce both (i) random numbers and (ii) inputs for the tables:

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

Variance Gamma Distribution

Definition

Variance gamma distribution is related to one of the most popular examples of Lévy processes, the Variance gamma process. There are several ways to its formulation, either as a difference of two gamma processes, which are, clearly, independent, or as a (geometric) Brownian motion subordinated by a gamma process, ie the gamma process gives an internal time for the overall process. Since VG distribution can be obtained by putting a random number from Gamma distribution into Gaussian distribution, we can denote it as VG[g(t,ν);θ,ϑ], ie the Brownian motion with mean θ - it drives the skewness, and standard deviation ϑ within a random time specified by Gamma distribution with unique mean and variance ν - it determines the kurtosis. Hence, defining the Brownian motion on the basis of standard normal distribution variates, ε, within standard clock time t, , we get VG process as follows: .

The characteristic function of VG distribution is (assuming the triple of parameters ν, θ, and ϑ):

It can be obtained by combining the characteristic functions of gamma and normal distribution. Similarly, when we put together the probability distribution function of Gamma and Normal distribution, or a distribution function:

we get a probability distribution function of VG distribution:

In order to precalculate particular moments we should utilize the characteristic function. That is, k-th moment of random variable X∈VG: . This can be easily transformed to centered moments and also to normalized moments (skewness and kurtosis).

We plot the probability distribution function bellow.

Note: with plain settings the first function is drawn with blue, the second with violet.

Algorithms

This part is under construction...

Since there is no predefined algorithm in Mathematica to obtain VG random numbers, we must define some, at elast on the basis of normally and gamma distributed variates. However, it is very slow.

PMC (Plain Monte Carlo)

In order to realize the simplest method, we need to utilize adjusted PMCrannum algorithm. Alternatively, we could put VGrandom algorithm, however, it requires more time.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

AVMC (Monte Carlo with antithetic variates)

The most intuitive utilization of antithetic variates approach to MC simulation is to create a pair of normally distributed variates, keeping the rest of the original PMC algorithm intact.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

LHSMC (Latin Hypercube Monte Carlo)

The most intuitive utilization of antithetic variates approach to MC simulation is to create a pair of normally distributed variates, keeping the rest of the original PMC algorithm intact.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

AV-LHSMC (Latin Hypercube Sampled Monte Carlo with antithetic variates)

The most intuitive utilization of antithetic variates approach to MC simulation is to create a pair of normally distributed variates, keeping the rest of the original PMC algorithm intact.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

SSMC (Stratified Monte Carlo on the basis of VG distribution function)

The most intuitive utilization of antithetic variates approach to MC simulation is to create a pair of normally distributed variates, keeping the rest of the original PMC algorithm intact.

Now we map the Histogram function to vectors of random numbers generated for particular k.

To draw the QQ plot we first specify self-defined function to evaluate the CDF for each random number from the data set, specifing the DataRange on the x-axes. Next, we sort the original series of data and map the QQplot function to them.

In order to complete the comparison, we can evaluate the Kolmogorov-Smirnov test, ie we look for the widest distance between the simulated distribution function and the teoretical one.

Inverse Gaussian Distribution

Generally, we proceed in the same way as for the Gamma distribution, though, prespecified algorithm InverseGaussianDistribution can be used.

Normal Inverse Gaussian Distribution

Any of the algorithms that were stated above for Variance Gamma distribution can be used if Gamma distribution is replaced by Inverse Gaussian Distribution.

Plain vanilla call

Note: In this section we provide a sample for estimation of option price by MC approach under BS model. Particular models and approach differes only when error is required to be estimated, too. Otherwise, it is sufficient to replace suitable modul (see above).

Payoff function

In[25]:=

In[33]:=

In[34]:=

Out[34]=

BS model

Inputs (Vstupy)

Here we define basic input data: spot price of the underlying asset, exercise price, time to maturity, riskless rate, dividend rate (costs of cary), volatility of undelying returns, length of one time interval:

In[37]:=

Within the Black and Scholes model (BS model) the underlying asset price follows the Geometric Brownian motion: .

This gives us:

Analytical pricing (Analytické ocenění)

The first formula is the Black and Scholes call option pricing formula, which works also for the maturity time (it includes the payoff function). The second and the third formulas calculate the arguments to the distribution function of the normal distribution, which is defined in the last formula.

In[35]:=

For a given set of input data, the option price is as follows.

In[38]:=

Out[38]=

MC - general module

The simplest way to estimate the option price is to apply any of the modules defined above and obtain a vector of N random numbers. After that, terminal price of the underlying asset is calculated and the payoff function is evaluated. It can be done as follows.

In[49]:=

Such a definition of the module allows us to estimate the option price (the mean) and the error of the estimation. The next line gives us the resutls for N = scenarios, including its error and total time costs, for i=2,3,4,5,6.

In[51]:=

Out[51]//MatrixForm=

Created with Wolfram Mathematica 7.0