I find options fascinating because they deal with the abstract ideas of volatility and correlation, both of which are unobservable and can often seem like wild animal spirits (take the current stock market as an example). Understanding these subtle concepts is never easy, but it is essential in pricing some of the more exotic options which involve multiple underlying stocks. To set the scene, let’s pretend that your neighbor wants to make a bet with you where he will pay you $100 if Google (GOOG) and Apple (APPL) are above 500 and 240 respectively after 1 year, but you have to pay him $25 today. How would we determine if $25 is a good deal or not?

The Model

The first things we need to do is to define the dynamics for all the tradable assets. We have at our disposal a money market account, as well as the stocks of Google and Apple (it is not a coincidence that I picked companies that don’t pay dividends). Below describes how we will model the dynamics of these assets.



In this model, is the price of a share in the money market account, and and are the stock prices. Notice that we have correlated the Brownian Motions such that the stock prices will tend to move up and down together just like in the real world. To make this model easier to work with, we will first convert the correlated Brownian Motions into independent Brownian Motions by defining the following.



By doing this, we can rewrite the ‘real world’ dynamics as

Calibration

Writing down a bunch of stochastic differential equations was pretty mechanical, but now comes the artsy part of modeling. The first thing we need is some stock price data, which can easily be downloaded from Yahoo Finance. This data gives us the stock price over discreet time intervals (I picked daily close prices). To use this data, lets define the discreet return as the solution to the following , in which case we can write . By assuming the drift and volatility terms are constants, we can use our model of the stock price dynamics to write the discreet returns as the following.



Where is the discreet time interval. By doing this, we can show that the discreet returns are distributed iid normal.



We can use the iid normal distribution to calibrate all the parameters. are simply the sample standard deviations of the discreet returns divided by , is the sample correlation of the discreet returns, and are the sample means of the discreet returns divided by (although I should note that we don’t actually need to know to price this option).

The Data

Now that we have the math worked out for calibrating the model, let’s take a look at some real data and see if the data agrees with our assumptions. First we need to read the data into R. The data I used can be found here. (I used R’s scan function since I’m lazy).

## Scan in close prices
GOOG <- scan()
AAPL <- scan()

Let’s take a look at the time series of stock prices.

## regular stock price chart
par(mfrow=c(2,1))
plot(GOOG, type="l", main="GOOG closing Price")
plot(AAPL, type="l", main="AAPL closing Price")

Now calculate and plot the daily returns.

## convert to returns
n <- length(GOOG)
GOOG.r <- log(GOOG[2:n]/GOOG[1:(n-1)])
AAPL.r <- log(AAPL[2:n]/AAPL[1:(n-1)])
 
## plot returns
matplot(cbind(GOOG.r, AAPL.r), type="l", lty=1, main="Daily Returns", ylab="Daily Return")
abline(h=0, col="blue")
legend("bottomleft", c("GOOG", "AAPL"), lty=1, col=c(1,2))

To check the iid normal assumption, we should do an autocorrelation plot and a QQ normal plot.

## autocorrelation plot
par(mfrow=c(2,2))
acf(GOOG.r)
acf(AAPL.r)
 
## QQ normal plots
qqnorm(GOOG.r, main="GOOG")
qqline(GOOG.r, col="red")
qqnorm(AAPL.r, main="AAPL")
qqline(AAPL.r, col="red")

Lastly, lets also take a look at the joint distribution of the returns.

## pairs plot
plot(GOOG.r, AAPL.r, main="Daily Returns")

The ACF plot suggests that the daily returns are iid, but they certainly have heavier tails than what would be produced by a normal distribution. In order to keep the model simple, we’ll assume normality of the returns holds, and give ourselves a ‘buffer’ when we calculate the fair price. Finally, time for the calibration!

GOOG vol: 0.2262006 AAPL vol: 0.2756421 Rho: 0.5413732

## calculate annualized vol
GOOG.vol <- sd(GOOG.r)/sqrt(1/252)
AAPL.vol <- sd(AAPL.r)/sqrt(1/252)
 
## calculate rho
rho <- cor(GOOG.r, AAPL.r)
 
cat("GOOG vol:", GOOG.vol, "n")
cat("AAPL vol:", AAPL.vol, "n")
cat("Rho:", rho, "n")

Price Using Simulation

We are now at the point where we can price this option. As usual, the price of the option is the expected discounted payoff under the risk neutral measure. This particular option can be priced as where the risk neutral dynamics of the stock prices are the following.



Where are independent Brownian Motions under the risk neutral measure. We will estimate the expectation through Monte Carlo simulation under a Euler discretization scheme. Rather than assume a constant interest rate, lets make things a tad more interesting by modeling the interest rate using the Vasicek Model (I’ll assume a calibration, but we could have calibrated it to zero coupon bond prices).



The following R code simulates the risk neutral dynamics of this model and estimates the expectation. By doing so, we find that the fair price of this option is $0.31 per $1 of notional. Thus, buying it from our neighbor for $25 seems like a deal if we think the $6 difference is a sufficient buffer to cover the simplifying assumptions we made.

trials <- 10000     # simulation trials
T <- 1              # time until expiration (years)
m <- 200            # subintervals
dt <- T/m           # time per subinterval (years)
 
## set interest rate model parameters
k <- 0.4
theta <- 0.05
beta <- 0.03
 
## set strikes
K1 <- 500
K2 <- 240
 
## create storage vectors
s1 <- rep(0,m+1)       # stock price path 1
s2 <- rep(0,m+1)       # stock price path 2
r <- rep(0,m+1)        # interest rate
c <- rep(0, trials)    # payoff
 
## set initial stock prices and interest rate
s1[1] <- 500
s2[1] <- 240
r[1]  <- 0.01
 
## begin simulation
for(j in 1:trials){
  for(i in 2:(m+1)){
    z1 <- rnorm(1,0,1)
    z2 <- rnorm(1,0,1)
    z3 <- rnorm(1,0,1)
 
    ds1 <- s1[i-1]*(r[i-1]*dt + GOOG.vol*sqrt(dt)*z1)
    ds2 <- s2[i-1]*(r[i-1]*dt + AAPL.vol*sqrt(dt)*(rho*z1 + sqrt(1-rho^2)*z2))  
    dr  <- k*(theta - r[i-1])*dt + beta*sqrt(dt)*z3
 
    s1[i] <- s1[i-1] +  ds1
    s2[i] <- s2[i-1] +  ds2
    r[i]  <- r[i-1]  +  dr
  }
  ss <- sum(r[2:(m+1)]*dt)
  c[j] <- ifelse(s1[m+1]>K1 && s2[m+1]>K2, exp(-ss), 0)
}
 
cat("Option Price Estimate:", round(mean(c),3), "n")
cat("Standard Error:", round(sd(c)/sqrt(trials),3), "n")