Steven R. Dunbar
Department of Mathematics
203 Avery Hall
University of Nebraska-Lincoln
Lincoln, NE 68588-0130
http://www.math.unl.edu
Voice: 402-472-3731
Fax: 402-472-8466

Stochastic Processes and
Advanced Mathematical Finance

__________________________________________________________________________

Monte Carlo Simulation of Option Prices

_______________________________________________________________________

Note: These pages are prepared with MathJax. MathJax is an open source JavaScript display engine for mathematics that works in all browsers. See http://mathjax.org for details on supported browsers, accessibility, copy-and-paste, and other features.

_______________________________________________________________________________________________

Rating

Rating

Mathematically Mature: may contain mathematics beyond calculus with proofs.

_______________________________________________________________________________________________

Section Starter Question

Section Starter Question

What is a 95% confidence interval in a statistical experiment?

_______________________________________________________________________________________________

Key Concepts

Key Concepts

  1. Monte Carlo methods (or Monte Carlo experiments) are mathematical solution methods that rely on repeated random sampling to obtain numerical results.
  2. Variance reduction techniques create a new Monte Carlo problem with the same answer as the original but with a smaller value for the variance to increase the efficiency of the Monte Carlo estimation. Many methods for variance reduction for Monte Carlo methods exist.
  3. The idea behind importance sampling is that some values of the simulation have more effect on the estimation than others. The probability distribution of the sample is carefully changed to give “important” outcomes more weight to achieve some goal such as variance reduction.

__________________________________________________________________________

Vocabulary

Vocabulary

  1. Monte Carlo methods (or Monte Carlo experiments) are a large and diverse class of solution methods using random sampling from a specified distribution to generate numerical results for a mathematical model.
  2. The 95% level Monte Carlo confidence interval for 𝔼 g(X) is
    n t0.975 s n,n + t0.975 s n.

  3. The quantity sn is called the standard error.
  4. Variance reduction techniques create a new Monte Carlo problem with the same answer as the original but with a smaller value for the variance to increase the efficiency of Monte Carlo estimation.
  5. An antithetic sample is a sample that uses an opposite value of f(x), being low when f(x) is high and vice versa, to reduce the variance.
  6. Importance sampling is another variance reduction technique. The idea behind importance sampling is that some values in a simulation have more influence on the estimation than others. The probability distribution is carefully changed to give “important” outcomes more weight.

__________________________________________________________________________

Mathematical Ideas

Mathematical Ideas

Confidence Intervals for Option Prices

Monte Carlo methods (or Monte Carlo experiments) are a large and diverse class of solution methods using random sampling from a specified distribution to generate numerical results for a mathematical model. Monte Carlo methods are a statistical solution technique, among the various mathematical techniques that move from a mathematical model to a relation among the quantities of interest.

With Geometric Brownian Motion as a stochastic process model for security prices, we can use that process for simulating option prices with a Monte Carlo method. We assume that a security price is modeled by Geometric Brownian Motion, so we know the lognormal probability distribution at the expiration time. Draw n pseudo-random numbers x1,,xn from the lognormal distribution modeling the stock price S. Then approximate a put option price as the present value of the expected value of the function g(S) = max(K S, 0), with the sample mean

𝔼 g(S) 1 n i=1ng(x i) = n.

The examples in this section feature put options rather than the call options illustrated previously for a reason that will become clear later when we modify the Monte Carlo simulation with importance sampling.

The Weak and Strong Laws of Large Numbers ensure us that this approximation is valid for large n, and improves as the number of samples increases. Even more is true, since the Central Limit Theorem implies that the sample mean n is approximately normally distributed with mean 𝔼 g(S) and variance Var g(S) n,

n N(𝔼 g(S) , Var g(S) n).

This says that the value n that we estimate with the Monte Carlo simulations will have a deviation from the true expected value 𝔼 g(S) that is on the order of the standard deviation divided by n. More precisely, recall that for the standard normal distribution |Z| < 1.96 0.95. Then by the shifting and rescaling properties of the normal distribution, we can construct a corresponding 95% confidence interval for the estimate n as

𝔼 g(S) 1.96 Var g(S) n , 𝔼 g(S) + 1.96 Var g(S) n .

We interpret the interval above as saying that the estimate n will be in this interval in about 95% of samples. That is, if 100 people perform the sampling experiment with different draws from the pseudo-random number generator, then about 95 of them will generate an expected value that falls in this interval, but about 5 of them will be unlucky enough to draw a sample that generates an expected value falling outside the interval.

A problem with obtaining the confidence interval is that the mean 𝔼 g(S) and the variance Var g(S) are both unknown. These are respectively estimated with the sample mean n and the sample variance

s2 = 1 n 1 i=1n(g(x i) n)2

to obtain the 95% level Monte Carlo confidence interval for 𝔼 g(S)

n tn1,0.975 s n,n + tn1,0.975 s n.

The quantity sn is called the standard error. Note that the standard error is dependent on the drawn numbers x1,x2,,xn so it is subject to sample variability too. In fact, the sample quantity

(n 𝔼 g(S)) sn

has a probability distribution known as Student’s t-distribution, so the 95% confidence interval limits of ±1.96 must be modified with the corresponding 95% confidence limits of the appropriate Student’s t-distribution.

Another problem with the Monte Carlo confidence interval is that the width of the interval decreases as 1n, which is not particularly fast. To decrease the width of the confidence interval by half requires 4 times as many sample variates.

Simplified Example

Now apply confidence interval estimation to the present value of a simple put option price for a simplified security. The simplified security has a stochastic differential equation dS = rSdt + σSdW with a starting price S(0) = z0 = 1 at t = 0 with risk-free interest rate r = σ22, and a standard deviation σ = 1. Then the stochastic process is S(t) = eW(t) where W(t) is standard Brownian Motion. (The parameters are financially unrealistic, but lead to the simple exponential.) For the simplified put option the strike price is K = 1, and the time to expiration is T t = T 0 = 1. Short programs to create 95% confidence intervals for 100, 1,000, and 10,000 samples are in the Scripts section. The random seed is set to the arbitrary value 2015 so that for replication of the Monte Carlo simulation the sample script always yields the same results. The true value comes from the Black-Scholes put option value formula in Put-Call Parity Principle.. Some typical results from the R script are in Table 1.


nn95% Monte Carlo confidence interval



100 0.1580783 (0.1206047, 0.1955519)
1000 0.1539137 (0.1425191, 0.1653084)
10000 0.1451896 (0.1416714, 0.1487078)

Table 1: Monte Carlo confidence intervals for the present value of a simple put option values. The present value of the put option from the Black-Scholes formula is 0.1446101.

Although the Monte Carlo simulation yields good results fairly easily, a disadvantage is that it applies for only one set of parameter values. To find values for a variety of parameters requires more computations, but even that does not demonstrate the functional relationship among the parameters and the values. The next chapter will derive a complete mathematical solution for the value of European call and put options. However, mathematical solution is not always possible for more complex security derivatives, and then Monte Carlo simulation is a good alternative.

Variance Reduction Techniques

Applying Monte Carlo estimation to a random variable with a large variance creates a confidence interval that is correspondingly large. The width of the confidence interval can be reduced by increasing the sample size, but the rate of reduction is n. Reducing the width of the confidence interval by half requires 4 times more sample points. Sometimes we can find a way to reduce s instead. To do this, we construct a new Monte Carlo problem with the same answer as our original problem but with a smaller value for s. Methods to do this are known as variance reduction techniques. Since Geometric Brownian Motion has a variance that grows exponentially, Monte Carlo estimates of quantities based on Geometric Brownian Motion are excellent candidates for variance reduction techniques.

Variance reduction techniques increase the efficiency of Monte Carlo estimation. Having methods to reduce variability with a given number of sample points, or for efficiency to achieve the same variability with fewer sample points, is desirable. Following are two methods that illustrate typical variance reduction techniques.

Antithetic Sampling

An antithetic sample is a sample that gives a opposite value of g(x), being low when g(x) is high and vice versa. Ordinarily we get an opposite g by sampling at a point that is somehow opposite to x. The idea of antithetic sampling can be applied when it is possible to find transformations of X that leave its probability distribution unchanged. For example, if X is normally distributed N(0, 1), then X̃ = X is normally distributed as well, that is the antithetic sample we will apply here. More generally, if X has a symmetric distribution about 0, then X̃ = X has the same distribution. Another antithetic sample occurs when X is uniformly distributed on [0, 1] so X̃ = 1 X has the same probability distribution.

The antithetic sampling estimate of a value μ = 𝔼 g(X) with a sample size n is

μ̂anti = 1 2n j=1n g(X j) + g(X̃j) .

Suppose the variance is σ2 = 𝔼 (g(X) μ)2 < . Then the variance in the antithetic sample of size n is

Var μ̂anti = Var 1 2n j=1ng(X j) + g(X̃j) = n 4n2 Var g(X) + g(X̃) = 1 4n Var g(X) + Var g(X̃) + 2 Cov g(X),g(X̃) = σ2 2n(1 + ρ).

Previously, we have often assumed different samples are independent so that the covariance term ρ is 0, reducing the number of terms in an expression. However, here is a situation where we deliberately create statistical dependence with negative covariance 1 ρ < 0 to reduce the variance in the antithetic sampling estimate.

Table 2 gives the results of antithetic sampling on the simplified put option example after a simple change of variables to make the sampling distribution normal. Table 2 shows that for the put option Monte Carlo estimation not only are the option value estimates much better with antithetic sampling, but the confidence intervals are reduced to less than half of the width compared with Table 1.


nn95%Monte Carlo confidence interval



100 0.1429556 (0.1267189, 0.1591923)
1000 0.1432171 (0.1383814, 0.1480529)
10000 0.1452026 (0.1437125, 0.1466926)

Table 2: Monte Carlo confidence intervals with antithetic sampling for the present value of a simple put option. The present value of the put option from the Black-Scholes formula is 0.1446101.

Note that if g(x) is a function with odd symmetry, for example g(x) = cx, and X is symmetric about 0, for example normally distributed with antithetic sample X̃ = X, then μ = 0. Furthermore, the statistical estimate μ̂anti would be always be 0 and in fact the estimate could be achieved with just one sample point and no variance. If g is not symmetric, as in option estimation with g(x) = max(K x, 0), we can still expect some improvement. The degree of improvement in antithetic sampling depends on the relative contribution of the odd and even parts of the function g, see [1] for details.

Importance Sampling

Importance sampling is another variance reduction technique. The idea behind importance sampling is that some values in a simulation have more influence on the estimation than others. The probability distribution is carefully changed to give “important” outcomes more weight. If “important” values are emphasized by sampling more frequently, then the estimator variance can be reduced. Hence, the key to importance sampling is to choose a new sampling distribution that “encourages” the important values.

Let f(x) be the density of the random variable, so we are trying to estimate

𝔼 g(x) =g(x)f(x)dx.

Estimate 𝔼 g(x) with respect to another strictly positive density h(x). Then easily

𝔼 g(x) =g(x)f(x) h(x) h(x)dx.

or equivalently, estimate

𝔼 g(Y )f(Y ) h(Y ) = 𝔼 g̃(Y )

where Y is a new random variable with density h(y) and g̃(y) = g(y)f(y)h(y). For variance reduction, we wish to determine a new density h(y) so that Var g̃(Y ) < Var g(X). Consider

Var g̃(Y ) = 𝔼 g̃(Y )2 𝔼 g̃(Y ) 2 =g2(x)f2(x) h(x) dx 𝔼 g(X) 2.

By inspection, we can see that we can make Var g̃(Y ) = 0 by choosing h(x) = g(x)f(x)𝔼 g(X). This is the ultimate variance reduction but it relies on knowing 𝔼 g(X), the quantity we are trying to estimate. It also requires that g(x) is strictly positive.

To choose a good importance sampling distribution requires some educated guessing. Importance sampling is equivalent to a change of variables, and like a change of variables often requires insight. Each instance of importance sampling depends on the function and the distribution. A specific example will illustrate importance sampling.

The example is again to calculate confidence intervals for the present value of a Monte Carlo estimate of a European put option price, where g(x) = max(K x, 0) and S is distributed as a geometric Brownian motion derived from dS = rSdt + σSdW with S(0) = z0. To set some simple parameters, we take the initial value z0 = 1, the risk free interest rate r = σ22, the standard deviation σ = 1, the strike price K = 1 and the time to expiration T = 1. The security process is then S(t) = eW(t). At T = 1 the process has a lognormal distribution with parameters μ = 0 and σT = 1. Thus we start with

0 max(K x, 0) 1 2πσxT exp(1 2(ln(x)σT)2)dx.

Using K = 1, σ = 1, and T = 1 after a first change of variable the integral we are estimating is

𝔼 g(S) = max(1 ex, 0)ex22 2π dx =0(1 ex)ex22 2π dx.

Introduce the exponential distribution ey2 with a second change of variables x = y. Then the expectation integral becomes

01 ey 2πy ey2 2 dy.

This change of probability distribution has the added advantage that now the sampling from the exponential distribution is concentrated where the function contributes most to the integral. This explains why the examples use put options, not call options.

Now apply this importance sampling to the confidence interval estimation for the present value of put option prices. Short programs to create the confidence intervals for 100, 1,000, and 10,000 samples are in the Scripts section. Some typical results from the R script are in Table 1.


nn95% Monte Carlo confidence interval



100 0.1465861 (0.1395778, 0.1535944)
1000 0.1449724 (0.1427713, 0.1471736)
10000 0.1443417 (0.1436435, 0.1436435)

Table 3: Monte Carlo confidence intervals with importance sampling for the present value of a simple put option. The present value of the put option from the Black-Scholes formula is 0.1446101.

Table 3 shows that the put option Monte Carlo estimates with importance sampling are closer to the Black-Scholes value than the simple Monte Carlo sampling. The confidence intervals are about one-fifth the width compared with Table 1.

Application to Data

Standard and Poor’s 500 Index, with symbol SPX, is a weighted index of 500 stocks. This index is designed to measure performance of the broad domestic economy through changes in the aggregate market value of 500 stocks representing all major industries. As such, it is a good candidate for a stock index modeled by Geometric Brownian Motion. The index started with a level of 10 for the 1941-43 base period.

To test Monte Carlo simulation on real data, use historical data on put options on the SPX that provides values for K, T t, and z0. Values for r and σ can be inferred from historical data.

The expectation to calculate for the Monte Carlo simulation of the value of a put option is

0 max(K x, 0) 1 2πσxt exp 1 2 ln(x) ln(z0) (r σ22)t σt 2 dx.

Then the expectation integral is equivalent to

0 max(K z 0e(rσ22)t+σtu, 0) + max(K z 0e(rσ22)tσtu, 0) 1 2π2u exp u 2 du.

This is now the importance sampling estimate of the value of a put option on the SPX index. The Scripts section has programs to evaluate confidence intervals for the put option value using real data.

Sources

This section is adapted from: Simulation and Inference for Stochastic Differential Equations, by Stephen Iacus, Sections 1.3 and 1.4, Springer-Verlag, 2008. The section on antithetic sample is adapted from Owen, Variance Reduction. and Lehna.. Interest rates are from Treasury.gov.

_______________________________________________________________________________________________

Algorithms, Scripts, Simulations

Algorithms, Scripts, Simulations

Algorithm

This script calculates the value of a put option on the Standard & Poor’s 500 stock index (SPX) in several ways. First the script uses the Black Scholes formula for a put option to calculate the theoretical value for a put option. Then the script uses Monte Carlo simulation based on the Geometric Brownian Motion with the SPX parameters to find an estimate for the put option value. Finally, the script uses Monte Carlo simulation with antithetic sampling and importance sampling variance reduction methods to refine the estimates. The ultimate purpose is to compare the various Monte Carlo methods with the put option value. The output will be a table of SPX put option value from the Black Scholes formula, Monte Carlo simulation with 95% confidence interval, antithetic sampling with 95% confidence interval and importance sampling with 95% confidence interval Set the total number of trials: n = 10000 Set the values of S, r, σ, K, T t Calculate the theoretical value of put option with Black-Scholes formula Fill the vector x[1..n] with random normal variates y1 = max(0,K S exp((r σ22)(T t) + σx(T t))) y2 = max(0,K S exp((r σ22)(T t) + σ(x)(T t))) The antithetic sample becomes y = (y1 + y2)2 Fill vector u[1..n] with random normal variates. The importance sample becomes u1 = (max(0,K S exp((r σ22)(T t) σT tx))+ max(0,K S exp((r σ22)(T t) + σT tx)))2πx. Apply the t-test to calculate mean and confidence intervals of the Monte Carlo simulations y1[1 : 100], y1[1 : 1000]and y1[1 : n]. Apply the t-test to calculate mean and confidence intervals of the antithetic samples y[1 : 100], y[1 : 1000] and y[1 : n] Apply the t-test to calculate mean and confidence intervals of the importance samples u1[1 : 100], u1[1 : 1000] and u1[1 : n]. Finally, collect and display a table of theoretical value and means and confidence intervals.

Scripts

R

R script for simplified option pricing with Monte Carlo simulation, antithetic sampling, and importance sampling..

1set.seed(2015) 
2n <- 10000 
3 
4S <- 1 
5sigma <- 1 
6r <- sigma^2/2 
7K <- 1 
8Tminust <- 1 
9 
10## the value of the put option by Black-Scholes 
11d2 <- (log(S/K) + (r - sigma^2/2) * (Tminust))/(sigma * sqrt(Tminust)) 
12d1 <- (log(S/K) + (r + sigma^2/2) * (Tminust))/(sigma * sqrt(Tminust)) 
13VP <- K * exp(-r * Tminust) * pnorm(-d2) - S * pnorm(-d1) 
14 
15x <- rlnorm(n)  #Note use of default meanlog=0, sdlog=1 
16y <- sapply(x, function(z) max(0, K - z)) 
17 
18mc100 <- t.test(exp(-r * Tminust) * y[1:100])  # first 100 simulations 
19mc1000 <- t.test(exp(-r * Tminust) * y[1:1000])  # first 1000 simulations 
20mcall <- t.test(exp(-r * Tminust) * y)  # all simulation results 
21 
22type <- c("Blacksholes", "Monte Carlo 100", "Monte Carlo 1000", "Monte Carlo all") 
23putestimate <- c(VP, mc100$estimate, mc1000$estimate, mcall$estimate) 
24putconfintleft <- c(NA, mc100$conf.int[1], mc1000$conf.int[1], mcall$conf.int[1]) 
25putconfintright <- c(NA, mc100$conf.int[2], mc1000$conf.int[2], mcall$conf.int[2]) 
26d <- data.frame(type, putestimate, putconfintleft, putconfintright) 
27 
28print(d) 
29 
30x <- rnorm(n) 
31y1 <- sapply(x, function(z) max(0, K - exp(z))) 
32y2 <- sapply(-x, function(z) max(0, K - exp(z))) 
33y <- (y1 + y2)/2 
34 
35anti100 <- t.test(exp(-r * Tminust) * y[1:100])  # first 100 simulations 
36anti1000 <- t.test(exp(-r * Tminust) * y[1:1000])  # first 1000 simulations 
37antiall <- t.test(exp(-r * Tminust) * y)  # all simulation results 
38 
39type <- c("Blacksholes", "Antithetic 100", "Antithetic 1000", "Antithetic all") 
40putestimate <- c(VP, anti100$estimate, anti1000$estimate, antiall$estimate) 
41putconfintleft <- c(NA, anti100$conf.int[1], anti1000$conf.int[1], antiall$conf.int[1]) 
42putconfintright <- c(NA, anti100$conf.int[2], anti1000$conf.int[2], antiall$conf.int[2]) 
43d <- data.frame(type, putestimate, putconfintleft, putconfintright) 
44 
45print(d) 
46 
47xN <- rnorm(n) 
48yN <- sapply(xN, function(z) max(0, K - exp(z))) 
49 
50mcN100 <- t.test(exp(-r * Tminust) * yN[1:100])  # first 100 simulations 
51mcN1000 <- t.test(exp(-r * Tminust) * yN[1:1000])  # first 1000 simulations 
52mcNall <- t.test(exp(-r * Tminust) * yN)  # all simulation results 
53 
54type <- c("Blacksholes", "Monte Carlo 100", "Monte Carlo 1000", "Monte Carlo all") 
55putestimate <- c(VP, mcN100$estimate, mcN1000$estimate, mcNall$estimate) 
56putconfintleft <- c(NA, mcN100$conf.int[1], mcN1000$conf.int[1], mcNall$conf.int[1]) 
57putconfintright <- c(NA, mcN100$conf.int[2], mcN1000$conf.int[2], mcNall$conf.int[2]) 
58d <- data.frame(type, putestimate, putconfintleft, putconfintright) 
59 
60print(d) 
61 
62xE <- rexp(n, rate = 1/2) 
63yE <- sapply(xE, function(z) max(0, K - exp(-sqrt(z)))/(sqrt(2 * pi) * sqrt(z))) 
64 
65mcE100 <- t.test(exp(-r * Tminust) * yE[1:100])  # first 100 simulations 
66mcE1000 <- t.test(exp(-r * Tminust) * yE[1:1000])  # first 1000 simulations 
67mcEall <- t.test(exp(-r * Tminust) * yE)  # all simulation results 
68 
69type <- c("Blacksholes", "Monte Carlo 100", "Monte Carlo 1000", "Monte Carlo all") 
70putestimate <- c(VP, mcE100$estimate, mcE1000$estimate, mcEall$estimate) 
71putconfintleft <- c(NA, mcE100$conf.int[1], mcE1000$conf.int[1], mcEall$conf.int[1]) 
72putconfintright <- c(NA, mcE100$conf.int[2], mcE1000$conf.int[2], mcEall$conf.int[2]) 
73d <- data.frame(type, putestimate, putconfintleft, putconfintright) 
74 
75print(d)

R script for SPX option pricing..

Octave

Octave script for SPX option pricing..

1n = 10000; 
2 
3S = 1614.96; # Standard and Poors 500 Index, on 07 / 01 / 2013 
4r = 0.008; # implied risk free interest rate, between 3 year and 5 year T - bill rate 
5sigma = 0.1827; # implied volatility 
6K = 1575; # strike price 
7Tminust = 110 / 365; # 07 / 01 / 2013 to 10 / 19 / 2013 
8 
9## the true value from the Black Scholes Formula for put option 
10numerd1 = log(S / K) + (r + sigma ^ 2 / 2) * (Tminust); 
11numerd2 = log(S / K) + (r - sigma ^ 2 / 2) * (Tminust); 
12d1 = numerd1 / (sigma * sqrt(Tminust)); 
13d2 = numerd2 / (sigma * sqrt(Tminust)); 
14part1 = S * (stdnormal_cdf(d1) - 1); 
15part2 = K * exp(- r * (Tminust)) * (stdnormal_cdf(d2) - 1); 
16VP = part1 - part2; 
17 
18x = stdnormal_rnd(n, 1); 
19y1 = max(0, K - S * exp((r - sigma ^ 2 / 2) * Tminust + sigma * x * sqrt(Tminust))); 
20y2 = max(0, K - S * exp((r - sigma ^ 2 / 2) * Tminust + sigma * (- x) * sqrt(Tminust))); 
21y = (y1 + y2) / 2; 
22 
23u = exprnd(1 / 0.5, [n, 1]); # note exponential parameter 
24tildeg = ... 
25  (... 
26  max(0, K - S * exp((r - sigma ^ 2 / 2) * Tminust - sigma * sqrt(Tminust) * sqrt(u))) + ... 
27  max(0, K - S * exp((r - sigma ^ 2 / 2) * Tminust + sigma * sqrt(Tminust) * sqrt(u)))... 
28  ) ./ sqrt(2 * pi * u); # note use of broadcasting 
29 
30function [m, confintleft, confintright] = ttest(data, confidence) 
31m = mean(data); 
32se = std(data) / sqrt(length(data)); 
33df = length(data) - 1; 
34q = tinv((1 + confidence) / 2, df); 
35confintleft = m - q * se; 
36confintright = m + q * se; 
37endfunction 
38 
39[mc100mu, mc100cil, mc100cir] = ttest(exp(- r * Tminust) * y1(1:100), 0.95); 
40[mc1000mu, mc1000cil, mc1000cir] = ttest(exp(- r * Tminust) * y1(1:1000), 0.95); 
41[mcallmu, mcallcil, mcallcir] = ttest(exp(- r * Tminust) * y1, 0.95); 
42 
43[anti100mu, anti100cil, anti100cir] = ttest(exp(- r * Tminust) * y(1:100), 0.95); 
44[anti1000mu, anti1000cil, anti1000cir] = ttest(exp(- r * Tminust) * y(1:1000), 0.95); 
45[antiallmu, antiallcil, antiallcir] = ttest(exp(- r * Tminust) * y, 0.95); 
46 
47[imp100mu, imp100cil, imp100cir] = ttest(exp(- r * Tminust) * tildeg(1:100), 0.95); 
48[imp1000mu, imp1000cil, imp1000cir] = ttest(exp(- r * Tminust) * tildeg(1:1000), 0.95); 
49[impallmu, impallcil, impallcir] = ttest(exp(- r * Tminust) * tildeg, 0.95); 
50 
51putestimate = [VP, ... 
52  mc100mu, mc1000mu, mcallmu, ... 
53  anti100mu, anti1000mu, antiallmu, ... 
54  imp100mu, imp1000mu, impallmu]; 
55putconfintleft = [NA, ... 
56  mc100cil, mc1000cil, mcallcil, ... 
57  anti100cil, anti1000cil, antiallcil, ... 
58  imp100cil, imp1000cil, impallcil]; 
59putconfintright = [NA, ... 
60  mc100cir, mc1000cir, mcallcir, ... 
61  anti100cir, anti1000cir, antiallcir, ... 
62  imp100cir, imp1000cir, impallcir]; 
63d = transpose([putestimate; putconfintleft; putconfintright])
Perl

Perl script for SPX option pricing..

1use PDL::NiceSlice; 
2use PDL::GSL::RNG; 
3$rng = PDL::GSL::RNG->new(taus); 
4$rng->set_seed( time() ); 
5use PDL::Constants qw(PI E); 
6use Statistics::Distributions; 
7 
8$n = 10000; 
9 
10$S = 1614.96;    # Standard and Poors 500 Index, on 07/01/2013 
11$r = 
12    0.008; # implied risk free interest rate, between 3 and 5 year T-bill rate 
13$sigma   = 0.1827;       # implied volatility 
14$K       = 1575;         # strike price 
15$Tminust = 110 / 365;    # 07/01/2013 to 10/19/2013 
16 
17sub pnorm { 
18    my ( $x, $sigma, $mu ) = @_; 
19    $sigma = 1 unless defined($sigma); 
20    $mu    = 0 unless defined($mu); 
21 
22    return 0.5 * ( 1 + erf( ( $x - $mu ) / ( sqrt(2) * $sigma ) ) ); 
23} 
24 
25$numerd1 = log( $S / $K ) + ( ( $r + $sigma**2 / 2 ) * ($Tminust) ); 
26$numerd2 = log( $S / $K ) + ( ( $r - $sigma**2 / 2 ) * ($Tminust) ); 
27$d1 = $numerd1 / ( $sigma * sqrt($Tminust) ); 
28$d2 = $numerd2 / ( $sigma * sqrt($Tminust) ); 
29$part1 = $S * ( pnorm($d1) - 1 ); 
30$part2 = $K * exp( -$r * $Tminust ) * ( pnorm($d2) - 1 ); 
31$VP = $part1 - $part2; 
32 
33$x = $rng->ran_gaussian( 1, $n );    # tested to here 
34$y1 = pdl( 
35    zeroes($n), 
36    $K - $S * exp( 
37        ( $r - $sigma**2 / 2 ) * $Tminust + $sigma * $x * sqrt($Tminust) 
38    ) 
39)->transpose->maximum; 
40$y2 = pdl( 
41    zeroes($n), 
42    $K - $S * exp( 
43        ( $r - $sigma**2 / 2 ) * $Tminust + $sigma * ( -$x ) * sqrt($Tminust) 
44    ) 
45)->transpose->maximum; 
46$y = ( $y1 + $y2 ) / 2; 
47 
48$u = $rng->ran_exponential( 1 / 0.5, $n );    # note exponential parameter 
49$tildeg1 = pdl( 
50    zeroes($n), 
51    $K - $S * exp( 
52        ( $r - $sigma**2 / 2 ) * $Tminust - $sigma * sqrt($Tminust) * sqrt($u) 
53    ) 
54)->transpose->maximum; 
55$tildeg2 = pdl( 
56    zeroes($n), 
57    $K - $S * exp( 
58        ( $r - $sigma**2 / 2 ) * $Tminust + $sigma * sqrt($Tminust) * sqrt($u) 
59    ) 
60)->transpose->maximum; 
61$tildeg = ( $tildeg1 + $tildeg2 ) / sqrt( 2 * PI * $u ); 
62 
63sub ttest { 
64    my ( $data, $confidence ) = @_; 
65    $confidence = 0.95 unless defined($confidence); 
66    ( $mean, $prms, $median, $min, $max, $adev, $rms ) = stats($data); 
67    $df = nelem($data) - 1; 
68    $q = Statistics::Distributions::tdistr( $df, ( 1 - $confidence ) / 2 ); 
69 
70# Note use of level (1-$confidence)/2) instead of percentile (1 + $confidence)/2 
71    $widthci      = $q * $rms / sqrt( nelem($data) ); 
72    $confintleft  = $mean - $widthci; 
73    $confintright = $mean + $widthci; 
74    return ( $mean, $confintleft, $confintright ); 
75} 
76 
77( $mc100mu, $mc100cil, $mc100cir ) = 
78    ttest( exp( -$r * $Tminust ) * $y1 ( 0 : 99 ), 0.95 ); 
79( $mc1000mu, $mc1000cil, $mc1000cir ) = 
80    ttest( exp( -$r * $Tminust ) * $y1 ( 0 : 999 ), 0.95 ); 
81( $mcallmu, $mcallcil, $mcallcir ) = 
82    ttest( exp( -$r * $Tminust ) * $y1, 0.95 ); 
83 
84( $anti100mu, $anti100cil, $anti100cir ) = 
85    ttest( exp( -$r * $Tminust ) * $y ( 0 : 99 ), 0.95 ); 
86( $anti1000mu, $anti1000cil, $anti1000cir ) = 
87    ttest( exp( -$r * $Tminust ) * $y ( 0 : 999 ), 0.95 ); 
88( $antiallmu, $antiallcil, $antiallcir ) = 
89    ttest( exp( -$r * $Tminust ) * $y, 0.95 ); 
90 
91( $imp100mu, $imp100cil, $imp100cir ) = 
92    ttest( exp( -$r * $Tminust ) * $tildeg ( 0 : 99 ), 0.95 ); 
93( $imp1000mu, $imp1000cil, $imp1000cir ) = 
94    ttest( exp( -$r * $Tminust ) * $tildeg ( 0 : 999 ), 0.95 ); 
95( $impallmu, $impallcil, $impallcir ) = 
96    ttest( exp( -$r * $Tminust ) * $tildeg, 0.95 ); 
97 
98$d = pdl( 
99    [ $VP,       "NaN",      "NaN" ], 
100    [ $mc100mu,  $mc100cil,  $mc100cir ], 
101    [ $mc1000mu, $mc1000cil, $mc1000cir ], 
102    [ $mcallmu,  $mcallcil,  $mcallcir ], 
103 
104    [ $anti100mu,  $anti100cil,  $anti100cir ], 
105    [ $anti1000mu, $anti1000cil, $anti1000cir ], 
106    [ $antiallmu,  $antiallcil,  $antiallcir ], 
107 
108    [ $imp100mu,  $imp100cil,  $imp100cir ], 
109    [ $imp1000mu, $imp1000cil, $imp1000cir ], 
110    [ $impallmu,  $impallcil,  $impallcir ] 
111); 
112 
113print $d;
SciPy

Python script for SPX option pricing..

1import scipy 
2from scipy.stats import norm 
3 
4n = 10000 
5 
6S = 1614.96  # Standard and Poors 500 Index, on 07/01/2013 
7r = 0.008  # implied risk free interest rate, between 3 year and 5 year T-bill rate 
8sigma = 0.1827  # implied volatility 
9K = 1575  # strike price 
10Tminust = 110. / 365.  # 07/01/2013 to 10/19/2013 
11 
12numerd1 = scipy.log(S / K) + (r + sigma ** 2 / 2) * Tminust 
13numerd2 = scipy.log(S / K) + (r - sigma ** 2 / 2) * Tminust 
14d1 = numerd1 / (sigma * scipy.sqrt(Tminust)) 
15d2 = numerd2 / (sigma * scipy.sqrt(Tminust)) 
16part1 = S * (norm.cdf(d1) - 1) 
17part2 = K * scipy.exp(-r * Tminust) * (norm.cdf(d2) - 1) 
18VP = part1 - part2 
19 
20x = norm.rvs(size=n) 
21y1 = scipy.maximum(0, K - S * scipy.exp((r - sigma ** 2 / 2) * Tminust 
22                   + sigma * x * scipy.sqrt(Tminust))) 
23y2 = scipy.maximum(0, K - S * scipy.exp((r - sigma ** 2 / 2) * Tminust 
24                   + sigma * -x * scipy.sqrt(Tminust))) 
25y = (y1 + y2) / 2 
26 
27from scipy.stats import expon 
28u = expon.rvs(size=n, scale=1. / 0.5) 
29tildeg = (scipy.maximum(0, K - S * scipy.exp((r - sigma ** 2 / 2) 
30          * Tminust - sigma * scipy.sqrt(Tminust) * scipy.sqrt(u))) 
31          + scipy.maximum(0, K - S * scipy.exp((r - sigma ** 2 / 2) 
32          * Tminust + sigma * scipy.sqrt(Tminust) * scipy.sqrt(u)))) \ 
33    / scipy.sqrt(2 * scipy.pi * u) 
34 
35 
36def ttest(data, confidence=0.95): 
37    m = scipy.mean(data) 
38    se = scipy.stats.sem(data) 
39    df = len(data) - 1 
40    h = se * scipy.stats.t.ppf((1 + confidence) / 2., df) 
41    return (m, m - h, m + h) 
42 
43 
44[mc100mu, mc100cil, mc100cir] = ttest(scipy.exp(-r * Tminust) 
45        * y1[0:99], 0.95) 
46[mc1000mu, mc1000cil, mc1000cir] = ttest(scipy.exp(-r * Tminust) 
47        * y1[0:999], 0.95) 
48[mcallmu, mcallcil, mcallcir] = ttest(scipy.exp(-r * Tminust) * y1, 
49        0.95) 
50 
51[anti100mu, anti100cil, anti100cir] = ttest(scipy.exp(-r * Tminust) 
52        * y[0:99], 0.95) 
53[anti1000mu, anti1000cil, anti1000cir] = ttest(scipy.exp(-r * Tminust) 
54        * y[0:999], 0.95) 
55[antiallmu, antiallcil, antiallcir] = ttest(scipy.exp(-r * Tminust) 
56        * y, 0.95) 
57 
58[imp100mu, imp100cil, imp100cir] = ttest(scipy.exp(-r * Tminust) 
59        * tildeg[0:99], 0.95) 
60[imp1000mu, imp1000cil, imp1000cir] = ttest(scipy.exp(-r * Tminust) 
61        * tildeg[0:999], 0.95) 
62[impallmu, impallcil, impallcir] = ttest(scipy.exp(-r * Tminust) 
63        * tildeg, 0.95) 
64 
65putestimate = [ 
66    VP, 
67    mc100mu, 
68    mc1000mu, 
69    mcallmu, 
70    anti100mu, 
71    anti1000mu, 
72    antiallmu, 
73    imp100mu, 
74    imp1000mu, 
75    impallmu, 
76    ] 
77putconfintleft = [ 
78    scipy.nan, 
79    mc100cil, 
80    mc1000cil, 
81    mcallcil, 
82    anti100cil, 
83    anti1000cil, 
84    antiallcil, 
85    imp100cil, 
86    imp1000cil, 
87    impallcil, 
88    ] 
89putconfintright = [ 
90    scipy.nan, 
91    mc100cir, 
92    mc1000cir, 
93    mcallcir, 
94    anti100cir, 
95    anti1000cir, 
96    antiallcir, 
97    imp100cir, 
98    imp1000cir, 
99    impallcir, 
100    ] 
101d = scipy.transpose(scipy.array([putestimate, putconfintleft, 
102                    putconfintright])) 
103 
104print(d)

__________________________________________________________________________

Problems to Work

Problems to Work for Understanding

  1. Experiment with various values of sample size in the Monte Carlo simulation of the simplified put option to determine how the width of the confidence interval depends on the sample size.
  2. Experiment with various random seeds and values of the parameter σ in the simplified put option model to find the relative errors of the basic Monte Carlo method, antithetic sampling and importance sampling.
  3. From Options. we know the value of a put option increases with increasing volatility σ. Additionally the standard deviation of the simplified security price increases exponentially with increasing σ. Experiment with various values of σ in the Monte Carlo simulation of the simplified put option to determine how the put option price and the width of the confidence intervals vary with sigma.
  4. Evaluate the following integrals numerically (that is, with numerical integration software, not Monte Carlo evaluation):
    1. 0 max(K x, 0) 1 2πσxT exp(1 2(ln(x)σT)2)dx;

    2. 0(1 ex)ex22 2π dx;

    3. 01 ey 2πy ey2 2 dy.

    Interpreting each as the expected value of a put option at expiration, find the present value i of each if the interest rate is r = σ22 with σ = 1 and T t = T 0 = 1.

  5. Change the scripts to provide descriptive output.
  6. Change the scripts to create a Monte Carlo simulation estimate of the value of a put option on the Standard and Poor’s 500 index (SPX) directly using 10, 000 random variates with a lognormal distribution. Use parameters S = 1614.96, r = 0.008, σ = 0.1827, K = 1575, T t = 110365. Also construct the 95% confidence interval for the estimate.
  7. Change the scripts to create a Monte Carlo simulation estimate of the value of a call option on the Standard and Poor’s 500 index (SPX) directly using 10, 000 random variates with a lognormal distribution. Use parameters S = 1614.96, r = 0.008, σ = 0.1827, K = 1575, T t = 110365.
  8. Provide a detailed explanation of the changes of variables that lead from
    0 max(K x, 0) 1 2πσxt exp 1 2 ln(x) ln(z0) (r σ22)t σt 2 dx.

    to

    0 max(K z 0e(rσ22)t+σtu, 0) + max(K z 0e(rσ22)tσtu, 0) 1 2π2u exp u 2 du.

__________________________________________________________________________

Books

Reading Suggestion:

References

[1]   Art Owen. Monte carlo theory, methods and examples. http://statweb.stanford.edu/ owen/mc/, 2013.

__________________________________________________________________________

Links

Outside Readings and Links:

__________________________________________________________________________

I check all the information on each page for correctness and typographical errors. Nevertheless, some errors may occur and I would be grateful if you would alert me to such errors. I make every reasonable effort to present current and accurate information for public use, however I do not guarantee the accuracy or timeliness of information on this website. Your use of the information from this website is strictly voluntary and at your risk.

I have checked the links to external sites for usefulness. Links to external websites are provided as a convenience. I do not endorse, control, monitor, or guarantee the information contained in any external website. I don’t guarantee that the links are active at all times. Use the links here with the same caution as you would all information on the Internet. This website reflects the thoughts, interests and opinions of its author. They do not explicitly represent official positions or policies of my employer.

Information on this website is subject to change without notice.

Steve Dunbar’s Home Page, http://www.math.unl.edu/~sdunbar1

Email to Steve Dunbar, sdunbar1 at unl dot edu

Last modified: Processed from LATEX source on August 8, 2016