Chapter 11. Monte Carlo Simulation and Options

In finance, we study the trade-off between risk and return. The common definition of risk is uncertainty. For example, when evaluating a potential profitable project, we have to predict many factors in the life of the project, such as the annual sales, price of the final product, prices of raw materials, salary increase of employees, inflation rate, cost of borrowing, cost of new equity, and economic status. For those cases, the Monte Carlo simulation could be used to simulate many possible future outcomes, events, and their various combinations. In this chapter, we focus on the applications of the Monte Carlo simulation to price various options.

In this chapter, we will cover the following topics:

• Generating random numbers from standard normal distribution and normal distribution

• Generating random numbers from a uniform distribution

• A simple application: estimate pi by the Monte Carlo simulation

• Generating random numbers from a Poisson distribution...

Generating random numbers from a standard normal distribution

Normal distributions play a central role in finance. A major reason is that many finance theories, such as option theory and applications, are based on the assumption that stock returns follow a normal distribution. It is quite often that we need to generate n random numbers from a standard normal distribution. For this purpose, we have the following two lines of code:

>>>import scipy as sp
>>>x=sp.random.standard_normal(size=10)

The basic random numbers in SciPy/NumPy are created by Mersenne Twister PRNG in the numpy.random function. The random numbers for distributions in numpy.random are in cython/pyrex and are pretty fast. To print the first few observations, we use the print() function as follows:

>>>print x[0:5]
[-0.55062594 -0.51338547 -0.04208367 -0.66432268  0.49461661]
>>>

Alternatively, we could use the following code:

>>>import scipy as sp
>>>x=sp.random.normal...

Generating random numbers from a uniform distribution

When we plan to randomly choose m stocks from n available stocks, we could draw a set of random numbers from a uniform distribution. To generate 10 random numbers between one and 100 from a uniform distribution, we have the following code. To guarantee that we generate the same set of random numbers, we use the seed() function as follows:

>>>import scipy as sp
>>>sp.random.seed(123345)
>>>x=sp.random.uniform(low=1,high=100,size=10)

Again, low, high, and size are the three keywords for the three input variables. The first one specifies the minimum, the second one specifies the high end, while the size gives the number of the random numbers we intend to generate. The first five numbers are shown as follows:

>>>print x[0:5]
[ 30.32749021  20.58006409   2.43703988  76.15661293  75.06929084]
>>>

Using simulation to estimate the pi value

It is a good exercise to estimate pi by the Monte Carlo simulation. Let's draw a square with 2R as its side. If we put the largest circle inside the square, its radius will be R. In other words, the areas for those two shapes have the following equations:

By dividing equation (4) by equation (5), we have the following result:

In other words, the value of pi will be 4* Scircle/Ssquare. When running the simulation, we generate n pairs of x and y from a uniform distribution with a range of zero and 0.5. Then we estimate a distance that is the square root of the summation of the squared x and y, that is,. Obviously, when d is less than 0.5 (value of R), it will fall into the circle. We can imagine throwing a dart that falls into the circle. The value of the pi will take the following form:

The following graph illustrates these random points within a circle and within a square:

The Python program to estimate the value of pi is presented as follows:

import...

Generating random numbers from a Poisson distribution

To investigate the impact of private information, Easley, Kiefer, O'Hara, and Paperman (1996) designed a (PIN) Probability of informed trading measure that is derived based on the daily number of buyer-initiated trades and the number of seller-initiated trades. The fundamental aspect of their model is to assume that order arrivals follow a Poisson distribution. The following code shows how to generate n random numbers from a Poisson distribution:

import scipy as sp
import matplotlib.pyplot as plt
x=sp.random.poisson(lam=1, size=100)
#plt.plot(x,'o')
a = 5. # shape
n = 1000
s = np.random.power(a, n)
count, bins, ignored = plt.hist(s, bins=30)
x = np.linspace(0, 1, 100)
y = a*x**(a-1.)
normed_y = n*np.diff(bins)[0]*y
plt.plot(x, normed_y)
plt.show()  

Selecting m stocks randomly from n given stocks

Based on the preceding program, we could easily choose 20 stocks from 500 available securities. This is an important step if we intend to investigate...

Bootstrapping with/without replacements

Assume that we have the historical data, such as price and return, for a stock. Obviously, we could estimate their mean, standard deviation, and other related statistics. What are their expected annual mean and risk next year? The simplest, maybe naïve way is to use the historical mean and standard deviation. A better way is to construct the distribution of annual return and risk. This means that we have to find a way to use historical data more effectively to predict the future. In such cases, we could apply the bootstrapping methodology. For example, for one stock, we have its last 20-year monthly returns, that is, 240 observations.

To estimate next year's 12 monthly returns, we need to construct a return distribution. First, we choose 12 returns randomly from the historical return set without replacements and estimate their mean and standard deviations. We repeat this procedure 5,000 times. The final output will be our return-standard distribution...

Distribution of annual returns

It is a good application to estimate annualized return distribution and represent it as a graph. To make our exercise more meaningful, we download Microsoft 's daily price data. Then, we estimate its daily returns and convert them into annual ones. Based on those annual returns, we generate its distribution by applying bootstrapping with replacements 5,000 times as shown in the following code:

from matplotlib.finance import quotes_historical_yahoo
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
# Step 1: input area
ticker='MSFT'          # input value 1
begdate=(1926,1,1)     # input value 2
enddate=(2013,12,31)   # input value 3
n_simulation=5000      # input value 4
# Step 2: retrieve price data and estimate log returns
x=quotes_historical_yahoo(ticker,begdate,enddate,asobject=True,adjusted=True)
logret = log(x.aclose[1:]/x.aclose[:-1])
# Step 3: estimate annual returns
date=[]
d0=x.date
for i in range(0,size(logret)):
    date.append...

Simulation of stock price movements

We mentioned in the previous sections that in finance, returns are assumed to follow a normal distribution, whereas prices follow a lognormal distribution. The stock price at time t+1 is a function of the stock price at t, mean, standard deviation, and the time interval as shown in the following formula:

In this formula, is the stock price at t+1, is the expected stock return, is the time interval (), T is the time (in years), n is the number of steps, ε is the distribution term with a zero mean, and σ is the volatility of the underlying stock. With a simple manipulation, equation (4) can lead to the following equation that we will use in our programs:

In a risk-neutral work, no investors require compensation for bearing risk. In other words, in such a world, the expected return on any security (investment) is the risk-free rate. Thus, in a risk-neutral world, the previous equation becomes the following equation:

If you want to learn more about the risk...

Finding an efficient portfolio and frontier

In this section, we show you how to use the Monte Carlo simulation to generate returns for a pair of stocks with known means, standard deviations, and correlation between them. By applying the maximize function, we minimize the portfolio risk of this two-stock portfolio. Then, we change the correlations between the two stocks to illustrate the impact of correlation on our efficient frontier. The last one is the most complex one since it constructs an efficient frontier based on n stocks.

Finding an efficient frontier based on two stocks

The following program aims at generating an efficient frontier based on two stocks with known means, standard deviations, and correlation. We have just six input values: two means, two standard deviations, the correlation (), and the number of simulations. To generate the correlated y1 and y2 time series, we generate the uncorrelated x1 and x2 series first. Then, we apply the following formulae:

Another important issue...

Geometric versus arithmetic mean

In the next section, we discuss long-term return forecasting. Since we apply the weighted arithmetic and geometric means, we need to familiarize ourselves with the geometric mean first. For n returns (,, ) their arithmetic and geometric means are defined as follows:

In this formula, Ri is the stock's ith return. For an arithmetic mean, we could use the mean() function. Most of the time, the arithmetic mean is used in our estimations because of its simplicity. Since geometric means consider the time of values, it is considered to be more accurate for returns' estimation based on historical data. One important feature is that the geometric mean is smaller than its corresponding arithmetic mean unless all input values, such as all returns, are all the same. Because of this feature, many argue that using arithmetic means to predict future returns would lead to an overestimation. In contrast, geometric means would lead to an underestimation. Since the geometric...

Long-term return forecasting

Many researchers and practitioners argue that a long-term return forecast would be overestimated if it is based on the arithmetic mean of the past returns and underestimated based on a geometric mean. Using 80 years' historical returns to forecast the next 25-year future return, Jacquier, Kane, and Marcus (2003) suggest the following weighted scheme:

The following program reflects equation (12):

from matplotlib.finance import quotes_historical_yahoo
import numpy as np
import pandas as pd
ticker='IBM'           # input value 1
begdate=(1926,1,1)     # input value 2
enddate=(2013,12,31)   # input value 3
n_forecast=15.         # input value 4

def geomean_ret(returns):
    product = 1
    for ret in returns:
        product *= (1+ret)
    return product ** (1.0/len(returns))-1

x=quotes_historical_yahoo(ticker,begdate,enddate,asobject=True,adjusted=True)
logret = log(x.aclose[1:]/x.aclose[:-1])
date=[]
d0=x.date
for i in range(0,size(logret)):
    date.append(d0...

Pricing a call using simulation

After knowing the terminal prices, we could estimate the payoff for a call if the exercise price is given. The mean of those discounted payoffs using the risk-free rate as our discount rate will be our call price. The following code helps us estimate the call price:

from scipy import zeros, sqrt, shape
import scipy as sp
S0 = 40.      # stock price at time zero
X=   40.      # exercise price
T =0.5        # years
r =0.05       # risk-free rate
sigma = 0.2   # volatility (annual)
n_steps=100.          # number of steps
sp.random.seed(12345) # fix those random numbers
n_simulation = 5000   # number of simulation
dt =T/n_steps
call = zeros([n_simulation], dtype=float)
x = range(0, int(n_steps), 1)
for j in range(0, n_simulation):
    sT=S0
for i in x[:-1]:
        e=sp.random.normal()
        sT*=exp((r-0.5*sigma*sigma)*dt+sigma*e*sqrt(dt))
        call[j]=max(sT-X,0)
call_price=mean(call)*exp(-r*T)
print 'call price = ', round(call_price,3)

The estimated call...

Exotic options

Up to now, we discussed European and American options in Chapter 9, The Black-Scholes-Merton Option Model, which are also called vanilla options. One of the characters is path independent. On the other hand, exotic options are more complex since they might have several triggers relating to the determination of their payoffs. An exotic option could include nonstandard underlying instrument developed for particular investors, banks, or firms. Exotic options usually are traded over-the-counter (OTC). For exotic options, we don't have closed-form solutions, such as the Black-Scholes-Merton model. Thus, we have to depend on other means to price them. The Monte Carlo simulation is one of the ways to price many exotic options. In the next several subsections, we show how to price Asian options, digit options, and barrier options.

Using the Monte Carlo simulation to price average options

European and American options are path-independent options. This means that an option's payoff depends...

Barrier in-and-out parity

If we buy an up-and-out European call and an up-and-in European call, then the following parity should hold good:

The logic is very simple—if the stock price reaches the barrier, then the first call is worthless and the second call will be activated. If the stock price never touches the barrier, the first call will remain active, while the second one is never activated. Either way, one of them is active. The following Python program illustrates such scenarios:

def up_call(s0,x,T,r,sigma,n_simulation,barrier):
import scipy as sp
    import p4f
    n_steps=100.    
    dt=T/n_steps
    inTotal=0
    outTotal=0
    for j in range(0, n_simulation):
        sT=s0
        inStatus=False
        outStatus=True
for i in range(0,int(n_steps)):
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt))
if sT>barrier:
                outStatus=False
                inStatus=True
                #print 'sT=',sT
        #print 'j='...

Pricing lookback options with floating strikes

The lookback options depend on the paths (history) travelled by the underlying security. Thus, they are called path-dependent exotic options as well. One of them is named floating strikes. The payoff function of a call when the exercise price is the minimum price achieved during the life of the option is given as follows:

The Python code for this lookback option is shown as follows:

def lookback_min_price_as_strike(s,T,r,sigma,n_simulation):
    n_steps=100.    
    dt=T/n_steps
    total=0
    for j in range(n_simulation):
        min_price=100000.   # a very big number 
        sT=s
        for i in range(int(n_steps)):
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt))
if sT<min_price:
                min_price=sT
            #print 'j=',j,'i=',i,'total=',total
        total+=p4f.bs_call(s,min_price,T,r,sigma)
    return total/n_simulation    

Remember that the previous function needs two...

Using the Sobol sequence to improve the efficiency

When applying the Monte Carlo simulation to solve various finance related problems, we need to generate a certain number of random numbers. When the accuracy is very high, we have to draw a huge amount of such random numbers. For example, when pricing options, we use very small interval or a large number of steps to increase the number of decimal places of our final option prices. Thus, the efficiency of our Monte Carlo simulation would be a vital issue in terms of computational time and costs. This is especially true if we have a thousand options to price. One way to increase the efficiency is to apply a correct or better algorithm, that is, optimize our code. Another way is to use some special types of random number generators, such as the Sobol sequence.

Sobol sequences belong to the so-called low-discrepancy sequences, which satisfy the properties of random numbers but are distributed more evenly. Thus, they are usually called quasi random...

Summary

In this chapter, we discussed several types of distributions: normal, standard normal, lognormal, and Poisson. Since the assumption that stocks follow a lognormal distribution and returns follow a normal distribution is the cornerstone for option theory, the Monte Carlo simulation is used to price European options. Under certain scenarios, Asian options might be more effective in terms of hedging. Exotic options are more complex than the vanilla options since the former have no closed-form solution, while the latter could be priced by the Black-Scholes-Merton option model. One way to price these exotic options is to use the Monte Carlo simulation. The Python programs to price an Asian option and lookback options are discussed in detail.

In the next chapter, we will discuss various volatility measures, such as our conventional standard deviation, Lower Partial Standard Deviation (LPSD). Using the standard deviation of returns as a risk measure is based on a critical assumption that...

Exercises

1. Download daily price from Yahoo! Finance for DELL. Estimate daily returns and convert them into monthly returns. Assume its monthly returns follow a normal distribution. Draw a graph with the mean and standard deviation from the previous monthly returns.

2. Debug the following program:

import scipy as sp
S0 = 9.15 ;T =1;n_steps=10;mu =0.15;sigma = 0.2
n_simulation =  10
dt =T/n_steps    
S = sp.zeros([n_steps], dtype=float)
x = range(0, int(n_steps), 1)
for j in range(0, n_simulation):
    S[0]= S0
for i in x[:-1]:
        e=sp.random.normal()
        S[i+1]=S[i,j]+S[i]*(mu-0.5*pow(sigma,2))*dt+sigma*S[i]*sp.sqrt(dt)*e;
    plot(x, S)
figtext(0.2,0.8,'S0='+str(S0)+',mu='+str(mu)+',sigma='+str(sigma))
figtext(0.2,0.76,'T='+str(T)+', steps='+str(int(n_steps)))
title('Stock price (number of simulations = %d ' % n_simulation +')')
xlabel('Total number of steps ='+str(n_steps)))
ylabel('stock price')
show()

3. Write a Python program to price an Asian average price put based on the...