Monte Carlo Simulation in R with focus on Option PricingOjasvin SoodBlockedUnblockFollowFollowingJun 25In this blog, I will cover the basics of Monte Carlo Simulation, Random Number Distributions and the algorithms to generate them.
Finally I will also cover an application of Monte Carlo Simulation in the field of Option Pricing.
The whole blog focuses on writing the codes in R, so that you can also implement your own applications of Monte Carlo Simulation in R.
What are Monte Carlo methods?Monte Carlo methods basically refer to class of algorithms which use Randomness to give an estimate.
Let’s take an example to show thisTo give a numerical estimate of this integral of a function using Monte Carlo methods, one can model this integral as E[f(U)] where U is uniform random number in [0,1].
Generate n uniform random variables between [0,1].
Let those be U₁,U₂,…Uₙ with function values f(U₁), f(U₂),…f(Uₙ) respectively.
Numerical Estimate of IntegralThus we got an estimate of the integral using n uniform random samples between [0,1].
This estimate would converge to “actual” α as n converges to ∞.
Now we have looked at how we got an estimate, let us look at how do you generate different random number distributions.
Let us also look at an example involving Derivative Pricing.
First I will be defining some important Financial Terminologies.
Stock Price S(t) will be a function of time t.
We will be considering an European Option and price it.
European Option: It is a Financial Instrument which gives you the right but does not obligate you to buy/sell assets at a price ( Strike Price (K)) which is fixed at a previous point of time, lets say time 0.
Thus, if we go forward with the contract and agree to buy/sell an asset at a fixed price K at time T.
Our profit will beProfit = S(T)-K / K-S(T) respectivelyFrom now on lets only focus on buying an asset which is also called Call European OptionSo, if (S(T)-K) <0, then it is not in our interest to go forward with this contract.
Hence we don’t exercise our option i.
e.
we don’t buy the asset.
Hence, your payoff is-The value of PayoffBut, the question is, what is in it for the other party in this contract as you may or may not buy the asset, so the other party which currently holds the asset actually goes through with the contract as you also pay it some premium/fees initially.
Therefore the most important thing now is how much should this premium be ?You pay a premium amount at the beginning which is say C.
From this payoff we calculate profit as, Profit= Payoff — CThus, when we talk about derivative pricing, we try to estimate this C.
I will not go into the details now and cover it when we talk later about Geometric Brownian Motion, but the Stock Price variation is calculated using a Stochastic Differential Equation and upon solving it we get,S(T) Equation as a function of THere Z is another random number known as a Standard Normal Random variable.
We look in future on how to generate these random numbers but for now let us assume they can be generated by us.
Now we know the profit at a time T and also the value of Stock Price at time T.
How do we calculate it for time t=0.
Again with out going to the details the Premium C at time 0 is represented as E [e⁻ʳᵗ ( max ( S(t) — K ), 0 ) ].
Thus this is similar to our E[f(U)] which we saw earlier.
Hence, we again take help of Monte Carlo Simulation to give an estimate of this Payoff at time t=0.
Algorithm to Estimate European Option PayoffHence we repeatedly sample Z to find Stock price values and finally calculate European Option Premium C.
Generating Random DistributionsNow the only missing thing in previous cases is how would one generate a Uniform random, Normal random distributions.
We therefore look to cover algorithms to generate such Uniform random distributions and also methods to transform these to other distributions such as Normal Distributions.
The numbers that we will be generating in R for these simulations are actually called pseudorandom numbers which are not truly random but pretty good at mimicking genuine randomness.
We generate them using deterministic algorithms such that they approximate the properties of sequence of random numbers.
Uniform DistributionUniform Distribution is a continuous distribution abbreviated as U(a,b).
Where the probability density function isP.
D.
F of Uniform DistributionThe most popular method for generating such numbers is the Linear Congruential Generators.
The LCG is generated using-Linear Congruential Generator AlgorithmThus if we choose a x₀ also called the seed.
We can generate m-1 unique Uniform Random Numbers between [0,1] .
Though a good estimator must follow certain conditions for non zero c.
The above code produces Uniform Random Numbers.
Now let us look at Exponential, Normal Distribution and how we can generate them if we assume we already have a Uniform Random Number u.
Exponential DistributionExponential Distribution can be represented using the probability density functionP.
D.
F of Exponential DistributionTo generate it using uniform distribution we look at a Inverse Transform MethodInverse Transform EquationThus using a Uniform Random Number u and taking the inverse of Distribution function of Exponential Distribution we can generate a exponential random number X.
The proof of correctness for Inverse Transform Method is given asProof of Inverse Transform MethodHence we use this method to generate Exponential Distribution for which F⁻¹(U) is given asExponential Distribution using Inverse Transform MethodHere θ is exponential parameter 1/α.
Normal DistributionNormal distribution, also known as the Gaussian distribution, is a probability distribution that is symmetric about the mean, showing that data near the mean are more frequent in occurrence than data far from the mean.
The probability density function which depends on two parameters- mean and variance isP.
D.
F of Normal DistributionWe will look at a very popular method for generating standard Normal Distribution known as Box Muller Method.
The method is based on the fact that if we take a bivariate standard normal distribution Z₁, Z₂ ~ N (0, 1) then they follow two properties thatUsing the fact (ii) we generate the radius of the circle as R=sqrt(-2 log ( U₁ )) and the angle at which the point ( Z₁, Z₂ ) are located in the circle as V=2 π ( U₂ ) using a pair of Uniform Random Numbers U₁, U₂.
Therefore the pair of Z₁, Z₂ will be generated asLets look at the code and the plots nowIn the above code runif(5000) in R is used to generate 5000 Uniform Random Numbers in range [0,1].
Those are then used to generate finally 10000 Standard Normal Random Variables.
We can further improve on this using Marsaglia-Bray algorithm which I won’t cover in this series.
You can have a look at it on the link attached.
Also I would encourage you to code Marsaglia-Bray algorithm and check how fast it is as compared to Box-Muller algorithm.
Brownian MotionNow that the concepts on basics of Monte Carlo Simulation and various Random Distributions have been introduced lets focus on using Monte Carlo methods to simulate paths for various Stochastic Processes.
Standard Brownian Motion on [0,T] is a Stochastic Process { W(t), 0≤t≤T } which satisfies some properties such as-i) W(0)=0ii) For any k and any 0≤ t₁ ≤ t₂ ≤ ….
≤ tₖ ≤ T, the increments between any two successive W(tᵢ)-W(tᵢ-₁) are independent.
iii) The difference W(t)- W(s) ~ N(0,t-s) for any 0≤ s < t ≤ TAs a consequence of i) and iii) W(t) ~ N(0,t).
On the other hand, Brownian Motion which is non-standard will have two parameters just like Normal Distribution known as drift and diffusion.
Using W(t) we therefore give a Stochastic Differential Equation for any Brownian MotionSample Paths GenerationsSolving the SDE presented above we can write the equation in terms of X(tᵢ), μ(s), σ(s)Solution of Stochastic Differential Equation for Brownian MotionHence let us look at the code to generate paths where I have assumed μ and σ to be constantThe above code used the formula defined for Brownian Motion to generate 10 paths for it and used rnorm(1) a function defined in R to generate Standard Normal Random Variable.
The above image represents the 10 paths the code generated for the Standard Brownian Motion.
Geometric Brownian MotionA Stochastic process S(t) is a Geometric Brownian Motion if the differencelog(S(t))-log(S(0)) is a Brownian Motion.
Fundamentally, what that means in comparison to property ii) as discussed in Brownian Motion is that now,are independent of each other rather than just the successive differences as in Brownian Motion.
Hence, the Stochastic Differential Equation that we use now to represent the GBM isS.
D.
E for Geometric Brownian MotionThe solution to which is the same as what equation was used earlier to estimate the Stock Price where I had shown the example of using Monte Carlo Simulation.
Solution of S.
D.
E for Geometric Brownian Motionwhere W(t) is replaced withand thats how by using Monte Carlo Simulation we could also simulate the path of a Stock Price or a Geometric Brownian Motion.
For such simulation we again would have to discretize the time line into some N points to generate Stock Price at all such points.
Let us take initial Stock Price to be 100The plot looksHence we can use all such paths to finally get a S(T) and calculate the premium of European Option and finally give the average estimate of the European Option Price.
Next we will look at another Option known as Asian Option which would need the intermediate Stock Price values as well to calculate the Option Price.
Asian Option PricingIn an Asian option, the payoff depends on the average price of the underlying asset over a certain period of time as opposed to standard options where it depends only on the price at maturity.
Average Price Calculation for Asian OptionPayoff now will bePayoff of Asian OptionOur interest in simulating paths of a Stock Price was actually in pricing these options specially as the payoff depended on the path of the underlying asset S.
Let us look at the code for Asian Option PricingConclusionI hope you all get a fair introduction to not only Monte Carlo methods but also the field of Financial Engineering (Option Pricing).
Now you should be familiar with Monte Carlo methods, Derivative Pricing (European and Asian Options), Random Number Distributions (Uniform, Exponential and Normal Distributions) , basics of programming in R, Geometric Brownian Motion and its path generation.
.