European Options Valuation

Examples by Prof. Jacob Escobar, EGADE Business School, Tec de Monterrey

# First, we load some packages for prettier graphs
library(ggplot2) # already installed in Google Colab
# install.packages("gridExtra") # this isn't included in Google Colab
library(gridExtra)

Total Value

\(Total.value = Intrinsic.value + Time.value\)

Intrinsic value: the maximum of zero and the payoff from the option if it were exercised immediately. * Call: \(max(S – K, 0)\) * Put: \(max(K – S, 0)\)

Time value: The excess of an option’s value over its intrinsic value

Moneyness

  • At-the-money (ATM): S = K
  • In-the-money (ITM): Option has intrinsic value > 0
  • Out-of-the-money (OOTM): Option has instrinsic value = 0
##### EUROPEAN OPTIONS PAYOFF #####
# Assuming a $50 strike price (k)
k <- 50
# We define a vector of potential terminal spot prices (st)
st <- 0:100
st
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100
# We calculate the potential call and put options payoffs
# max() vs pmax()
c_po <- pmax(st - k, 0)
c_po
##   [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##  [51]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
##  [76] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## [101] 50
p_po <- pmax(___, ___)
p_po
po <- data.frame(st, c_po, p_po)
print(po[po$st >= 40 & po$st <= 60,], row.names = F)
##  st c_po p_po
##  40    0   10
##  41    0    9
##  42    0    8
##  43    0    7
##  44    0    6
##  45    0    5
##  46    0    4
##  47    0    3
##  48    0    2
##  49    0    1
##  50    0    0
##  51    1    0
##  52    2    0
##  53    3    0
##  54    4    0
##  55    5    0
##  56    6    0
##  57    7    0
##  58    8    0
##  59    9    0
##  60   10    0

Payoffs

Call * Long: \(+max(S - K, 0)\) * Short: \(-max(S - K, 0)\)

Put * Long: \(+max(K - S, 0)\) * Short: \(-max(K - S, 0)\)

# We graph the payoff of each combination of option type (call or put)
# and trade position (long or short)
# Notice that the payoffs of the short positions are the negative of the 
# payoffs of the long positions
lcpo <- ggplot(po, aes(x = st, y = c_po)) + 
  geom_line(color = "blue", size = 1) +
  ggtitle("Long Call Payoff")

lppo <- ggplot(po, aes(x = st, y = ___)) + 
  geom_line(color = "blue", size = 1) +
  ggtitle("Long Put Payoff")

scpo <- ggplot(po, aes(x = st, y = ___)) + 
  geom_line(color = "magenta", size = 1) + 
  ggtitle("Short Call Payoff")

sppo <- ggplot(po, aes(x = st, y = ___)) + 
  geom_line(color = "magenta", size = 1) +
  ggtitle("Short Put Payoff")

grid.arrange(lcpo, lppo, scpo, sppo, nrow = 2)

Net Profits

Call * Long: \(max(S - K, 0) - c\) * Short: \(c - max(S - K, 0)\)

Put * Long: \(max(K - S, 0) - p\) * Short: \(p -max(K - S, 0)\)

##### EUROPEAN OPTIONS NET PROFIT #####
# If we assume call and put costs, we can chart the net profits
call = 10
put = 10
c_np <- c_po - call
p_np <- p_po - put
np <- data.frame(st, c_np, p_np)
print(np[np$st >= 35 & np$st <= 65,], row.names = F)
##  st c_np p_np
##  35  -10    5
##  36  -10    4
##  37  -10    3
##  38  -10    2
##  39  -10    1
##  40  -10    0
##  41  -10   -1
##  42  -10   -2
##  43  -10   -3
##  44  -10   -4
##  45  -10   -5
##  46  -10   -6
##  47  -10   -7
##  48  -10   -8
##  49  -10   -9
##  50  -10  -10
##  51   -9  -10
##  52   -8  -10
##  53   -7  -10
##  54   -6  -10
##  55   -5  -10
##  56   -4  -10
##  57   -3  -10
##  58   -2  -10
##  59   -1  -10
##  60    0  -10
##  61    1  -10
##  62    2  -10
##  63    3  -10
##  64    4  -10
##  65    5  -10
lcnp <- ggplot(np, aes(x = st, y = c_po)) +
  geom_line(color = "blue", size = 1) + 
  geom_line(y = c_np, linetype = "dashed") +
  ggtitle("Long Call Payoff & Net Profit") +
  ylim(-10, 50)

lpnp <- ggplot(np, aes(x = st, y = p_po)) +
  geom_line(color = "blue", size = 1) +
  geom_line(y = p_np, linetype = "dashed") +
  ggtitle("Long Put Payoff & Net Profit") +
  ylim(-10, 50)

scnp <- ggplot(np, aes(x = st, y = -c_po)) +
  geom_line(color = "magenta", size = 1) +
  geom_line(y = -c_np, linetype = "dashed") +
  ggtitle("Short Call Payoff & Net Profit") +
  ylim(-50, 10)

spnp <- ggplot(np, aes(x = st, y = -p_po)) +
  geom_line(color = "magenta", size = 1)  +
  geom_line(y = -p_np, linetype = "dashed") +
  ggtitle("Short Put Payoff & Net Profit") +
  ylim(-50, 10)
  
grid.arrange(lcnp, lpnp, scnp, spnp, nrow = 2)

Options Valuation

How many Cash Flows can we receive from an option? How could we valuate the option using the Discounted Cash Flows (DCF) model? Would it be easy?

Gordon’s Dividend Model for Stock’s Valuation


Zero dividend growth: \[ P_0 = \frac{D_0}{k_e} \] Constant dividend growth: \[ P_0 = \frac{D_0(1+g)}{k_e - g} \]

Capturing Risk in the DCF


Risk in the denominator: \[P_0 = \frac{CF}{(1+r_{cf})}\]

Risk in the numerator: \[P_0 = \frac{E(CF)}{(1+r_{f})}\]


E(CF) is the expected Cash Flow or “Certainty Equivalent Cash Flow”.

Monte Carlo Simulation

  • Involves generating multiple (\(10^6\) or more) random values corresponding to the independent variable, based on an assumed probability distribution.
  • For stock options valuation we want to model the price of the underlying asset, so the random values are plugged in the corresponding valuation function.
  • We then average the potential values of all scenarios, and bring that expected value to the present.

Theoretical Distribution of Stock Prices

Given by:

\[ S_T = S_0e^{[(\mu - \frac{\sigma^2}{2})T + \sigma\epsilon\sqrt{T}]} \]


* \(\mu\) is the risk free rate (the “drift”) * \(\epsilon\) is a random number with Standard Normal Distribution ~\(N(0,1)\) (the “stochastic” component, Wiener process, Brownian motion, etc) * Given that we are assuming that Returns follow a Normal Distribution, then Prices follow a Lognormal Distribution (they are always positive).

Step-by-Step Computation

Input data

  • type: european call or put
  • s: underlying asset’s spot price
  • k: option’s strike price
  • t: option’s time to maturity in years (e.g. 3 months = 3/12 years)
  • v: underlying asset’s returns annualized volatility
  • r: annualized risk-free rate
  • m: number of simulations
# Input data
m = 100000
type = "c"
s = 52
k = 50
t = 0.50
v = 0.20
r = 0.10
# Generation of stochastic component
epsilon <- rnorm(m, 0, 1)
hist(epsilon)

# Price simulation
st <- s*exp((r - (v^2)/2)*___ + v*___*sqrt(___))
print(paste("Min St:", round(min(st), 2), 
            "/ Avg St:", round(mean(st), 2),
            "/ Max St:", round(max(st), 2)))
hist(st)
# Payoff estimations for call or put  
if(type == "___") po <- pmax(st - k, 0)
if(type == "p") po <- pmax(___, ___)
hist(po)
# Risk-neutral stock price at T (Future)
round(s*exp(r*t), 2)
## [1] 54.67
# Scenarios
scenarios <- data.frame(epsilon, st, po)
scenarios[sample(nrow(scenarios), 15),]
##           epsilon       st         po
## 33087 -0.21077178 52.53271  2.5327133
## 762   -0.43097774 50.92196  0.9219599
## 63416 -0.18641739 52.71396  2.7139601
## 95450 -0.18016200 52.76061  2.7606140
## 69384  0.30011612 56.46870  6.4687026
## 20577 -1.37555350 44.55441  0.0000000
## 29238 -0.07465799 53.55373  3.5537319
## 93991 -0.57600058 49.88822  0.0000000
## 43217  0.03129158 54.36220  4.3621976
## 95546  0.89886831 61.45855 11.4585457
## 84622 -1.83803270 41.73360  0.0000000
## 85970 -0.77544900 48.50072  0.0000000
## 3268   0.38587405 57.15772  7.1577249
## 87725  0.74927639 60.17202 10.1720165
## 79676 -1.01661004 46.87448  0.0000000
# Option price
op <- mean(___) * exp(___)
op

Custom Function

We can define a custom function so that we can apply the method to a variety of options quickly.

For ease of use, I include default values for all inputs in this didactic example, but in a real setting I would only include a default value for the number of simulations m.

It wouldn’t make sense to include default values for the other inputs since in those cases I prefer to get an error message than to make a wrong assumption.

# Function definition
mc <- function(m = 100000, type = "c", s = 52, k = 50,
               t = 0.50, v = 0.20, r = 0.10) {
                 
  epsilon <- rnorm(m, 0, 1)
  st <- s*exp((r - (v^2)/2)*t + v*epsilon*sqrt(t))
  
  if(type == "c") po <- pmax(st - k, 0)
  if(type == "p") po <- pmax(k - st, 0)
  
  op <- mean(po) * exp(-r*t)
  return(op)
}
# To verify it was created successfully
mc
## function(m = 100000, type = "c", s = 52, k = 50,
##                t = 0.50, v = 0.20, r = 0.10) {
##                  
##   epsilon <- rnorm(m, 0, 1)
##   st <- s*exp((r - (v^2)/2)*t + v*epsilon*sqrt(t))
##   
##   if(type == "c") po <- pmax(st - k, 0)
##   if(type == "p") po <- pmax(k - st, 0)
##   
##   op <- mean(po) * exp(-r*t)
##   return(op)
## }
# Function use with default values
mc()
## [1] 5.561245
# Test the function with your own values, for example:
mc(type = "p")
## [1] 1.131118

Black-Scholes

Given that European Options can only be exercised at maturity, the only thing that matters is the prices distribution at time \(T\).

This option is NOT path-dependant, that is, it only matters the final \(S_T\), not how it got there.

Hence, we can “skip” the simulation part by using Black-Scholes equations.

Equations



Call price:

\[ c = S_0 N(d_1) - Ke^{-rT} N(d_2) \]

Put price:

\[ p = Ke^{-rT} N(-d_2) - S_0 N(-d_1) \]

Where \(N(x)\) is the Standard Normal Cumulative Distribution Function (cdf) for x:

\[ d_1 = \frac{ln(\frac{S_0}{K})+(r+\frac{\sigma^2}{2}) T}{\sigma \sqrt{T}} \]

\[ d_2 = \frac{ln(\frac{S_0}{K})+(r-\frac{\sigma^2}{2}) T}{\sigma \sqrt{T}} \]

Shorcuts:

\[ d_2 = d_1 - \sigma \sqrt{T} \] \[ N(-d_1) = 1 - N(d_1) \] \[ N(-d_2) = 1 - N(d_2) \]

Step-by-Step Computation

Input Data

  • type: call or put
  • s: underlying asset’s spot price
  • k: option’s strike price
  • t: option’s time to maturity in years (e.g. 3 months = 3/12 years)
  • v: underlying asset’s returns annualized volatility
  • r: annualized risk-free rate
# Input data
type <- "c"
s <- 52
k <- 50
t <- 0.50
v <- 0.20
r <- 0.10

# Black-Scholes equations
d1 <- (log(___)+(r+(v^2)/2)*t)/(v*sqrt(t))
d2 <- ___

if(type == ___) {
  op <- s*pnorm(d1) - k*exp(-r*t)*pnorm(d2)
}
if(type == "p") {
  op <- k*exp(___)*pnorm(-d2) - s*pnorm(___)
}
op

Custom Function

We can define a custom function so that we can apply the method to a variety of options quickly.

For ease of use, I include default values for all inputs in this didactic example, but in a real setting it wouldn’t make sense. I prefer to get an error message than to make a wrong assumption.

bs <- function(type = "c", s = 52, k = 50, 
               t = 0.5, v = 0.20, r = 0.10) {

  op <- NA
  d1 <- (log(s/k)+(r+(v^2)/2)*t)/(v*sqrt(t))
  d2 <- d1 - v*sqrt(t)
  if(type == "c")   {
    op <- s*pnorm(d1) - k*exp(-r*t)*pnorm(d2)
  }
  if(type == "p") {
    op <- k*exp(-r*t)*pnorm(-d2) - s*pnorm(-d1)
  }
  return(___)
}
# To verify it was created successfully
bs
## function(type = "c", s = 52, k = 50, 
##                t = 0.5, v = 0.20, r = 0.10) {
## 
##   op <- NA
##   d1 <- (log(s/k)+(r+(v^2)/2)*t)/(v*sqrt(t))
##   d2 <- d1 - v*sqrt(t)
##   if(type == "c")    {
##     op <- s*pnorm(d1) - k*exp(-r*t)*pnorm(d2)
##   }
##   if(type == "p") {
##     op <- k*exp(-r*t)*pnorm(-d2) - s*pnorm(-d1)
##   }
##   return(op)
## }
# Function use with default values
bs()
## [1] 5.564811
# Test the function with your own values, for example:
bs(type = "p")
## [1] 1.126282

Extension for Income Paying Assets



Call price:

\[ c = S_0 e^{-qT} N(d_1) - Ke^{-rT} N(d_2) \]

Put price:

\[ p = Ke^{-rT} N(-d_2) - S_0 e^{-qT} N(-d_1) \]

Where \(N(x)\) is the Standard Normal Cumulative Distribution Function (cdf) for x:

\[ d_1 = \frac{ln(\frac{S_0}{K})+(r-q+\frac{\sigma^2}{2}) T}{\sigma \sqrt{T}} \]

\[ d_2 = d_1 - \sigma \sqrt{T} \]

Sensitivity Analysis

We can leverage the use of the custom functions we made, to perform quick and easy sensitivity analyses.

  1. Construct a vector for a range of values for a particular input (for example v or t).
  2. Provide that vector to the bs (or mc) function, to get a range of resulting option prices.

Expectations

Variable c_eur p_eur c_ame p_ame
\(S_0\) + - + -
\(K\) - + - +
\(T\) ? ? + +
\(\sigma\) + + + +
\(r\) + - + -
\(D\) - + - +

Volatility

# Construction of Volatility Vector
vv <- seq(from=0.05, to=0.50, by=0.05)
vv
##  [1] 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
# Call and Put Prices Vector: Using a For Loop
cv <- 0
pv <- 0
for (i in 1:length(vv)) {
  cv[i] <- bs(v = vv[i])
  pv[i] <- bs(type = "p",v = vv[i])
}
cv
##  [1] 4.441803 4.612016 5.027362 5.564811 6.161241 6.789106 7.435065 8.091937
##  [9] 8.755498 9.423072
pv
##  [1] 0.003274707 0.173486927 0.588833682 1.126281845 1.722712587 2.350576843
##  [7] 2.996535798 3.653408017 4.316968887 4.984543547
# Call and Put Values: Supplying the Volatility Vector to the Valuation Function
# This is the same as the above, but in a more compact form
cv <- bs(v = ___)
cv
pv <- bs(type = ___, v = ___)
pv
# Sensitivity Plot
ggplot(data.frame(vv, cv, pv), aes(x = vv)) +
  geom_line(aes(y = cv, color = "Call")) +
  geom_line(aes(y = pv, color = "Put")) +
  scale_color_manual(values=c("red", "blue")) +
  labs(title = "Option Price Sensitivity to Volatility", 
       color = "", y = "Option Price", x = "Volatility")

Time To Maturity

# We take the "function with vector input" approach:
tv <- seq(from = 0, to = 4, by = 0.25)
cv <- bs(t = ___)
cv
pv <- bs(type = ___, t = ___)
pv
# Sensitivity Plot
ggplot(data.frame(tv, cv, pv), aes(x = tv)) +
  geom_line(aes(y = cv, color = "Call")) +
  geom_line(aes(y = pv, color = "Put")) +
  scale_color_manual(values=c("red", "blue")) +
  labs(title = "Option Price Sensitivity to Time to Maturity", 
       color = "", y = "Option Price", x = "Time to Maturity")

Currently holding back the potential gains of the Put: * Call is ITM and Put is OOTM * Risk-free rate is very high in this example (10%) * If s = k and \r = 0, then c = p

Asian Options Valuation

These Asian Options ARE path-dependent options, contrary to European Options.

Asian Fixed Strike (K):

  • \(call.payoff = max(S_{avg(0,T)} – K, 0)\)
  • \(put.payoff = max(K – S_{avg(0,T)}, 0)\)

Asian Floating Strike (K):

  • \(call.payoff = max(S_T – S_{avg(0,T)}, 0)\)
  • \(put.payoff = max(S_{avg(0,T)} – S_T, 0)\)

In any case, the average payoff at \(T\) of all scenarios should be discounted to present value:

  • \(option.price = E(payoff)e^{(-rT)}\)

Theoretical Distribution of Stock Prices

Given by:

\[ S_t = S_{t-1}e^{[(\mu - \frac{\sigma^2}{2})\Delta t + \sigma\epsilon\sqrt{\Delta t}]} \]


* Notice that we are no longer simulating each scenario’s underlier’s price in a single step from \(S_0\) to \(S_{m,T}\), but in \(n\) smaller steps * Given that we are splitting \(T\) in \(n\) steps, then: \(\Delta t = T/n\) * This will still result in a matrix of \(m\) scenarios (\(m\) rows), and \(n\) steps (\(n+1\) columns, since “today’s price” is column 1 and we want to model \(n\) steps after today) * Once we get the \(m\) potential \(S_T\)s, we estimate the expected (average) payoff and discount it to the present using \(T\)

Step-by-Step Computation

Input Data

  • type: european call or put
  • s: underlying asset’s spot price
  • k: option’s strike price
  • t: option’s time to maturity in years (e.g. 3 months = 3/12 years)
  • v: underlying asset’s returns annualized volatility
  • r: annualized risk-free rate
  • m: number of simulations
  • n: number of steps
m <- 10
n <- 10
s <- 52
k <- 50
t <- 0.50
v <- 0.20
r <- 0.10
# Delta t
dt <- ___
dt
st <- matrix(s, ___, ___)
st
# Price-paths simulation
for (j in 2:(n+1)) {
  epsilon <- rnorm(___)
  st[, j] <- st[, ___] * exp((r - (v^2)/2)*___ + v*___*sqrt(___))  
}
st
# Average price per path (by rows)
stavg <- apply(st, 1, mean)
stavg
##  [1] 51.95408 52.08657 57.44872 45.59164 61.73270 53.02723 61.56224 53.50524
##  [9] 45.50706 56.92144
# apply is more generalizable, but for means and sums we have specific functions
# like: rowMeans( ), colMeans( ), rowSums( ) & colSums( )
# this does the same as the above:
stavg <- rowMeans(st)
stavg
##  [1] 51.95408 52.08657 57.44872 45.59164 61.73270 53.02723 61.56224 53.50524
##  [9] 45.50706 56.92144

Ahora calcularemos los precios de las opciones para Calls y Puts Europeas, Asiáticas Fixed Strike y Asiáticas Floating Strike, aprovechando que la parte de la simulación de precios es la misma.

# empty matrix
op <- matrix(NA,2,3)
rownames(op) <- c("call","put")
colnames(op) <- c("eur","afix","afloat")
op
##      eur afix afloat
## call  NA   NA     NA
## put   NA   NA     NA
# populating the matrix with the result for each model
op[1, 1] <- mean(pmax(st[, ___] - k, 0)) # european call
op[2, 1] <- mean(pmax(k - st[n+1], 0)) # european put
op[1, 2] <- mean(pmax(stavg - k, 0)) # asian fixed strike call
op[2, 2] <- mean(pmax(k - stavg, 0)) # asian fixed strike put
op[1, 3] <- mean(pmax(st[, n+1] - stavg, 0)) # asian floating strike call
op[2, 3] <- mean(pmax(stavg - st[, n+1], 0)) # asian floating strike put

# these are the option prices at time T
op
# discounting to present value (in a single step): THESE ARE THE OPTION PRICES
op <- op*exp(___)
op
# sorting the price-paths according to St (the last price achieved)
st <- st[order(st[, n+1]),]
st
##       [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
##  [1,]   52 51.87232 47.24279 48.33314 46.80882 46.44897 43.99790 43.18239
##  [2,]   52 51.94665 50.02070 47.72998 45.65964 44.87377 43.31300 40.98308
##  [3,]   52 52.96943 54.20590 57.10936 55.90093 51.42392 50.96957 51.02916
##  [4,]   52 52.87624 55.65426 51.14424 53.75084 53.41495 53.75274 52.39285
##  [5,]   52 52.27435 53.41634 54.45457 56.28373 57.27959 57.48364 53.98322
##  [6,]   52 53.38779 56.34901 54.38466 56.99528 62.15878 63.01292 61.21479
##  [7,]   52 54.24976 50.93970 51.96717 52.42756 52.10473 53.35738 52.70991
##  [8,]   52 56.66804 54.78750 55.41555 58.69628 63.18116 65.57286 68.70031
##  [9,]   52 51.54630 54.90818 60.12866 63.45571 63.08581 65.88367 67.70385
## [10,]   52 52.02219 51.89943 53.74044 51.80452 54.52734 58.93521 59.86333
##           [,9]    [,10]    [,11]
##  [1,] 40.62498 40.37860 40.61808
##  [2,] 40.88924 41.45263 41.70894
##  [3,] 48.74215 49.31186 47.83265
##  [4,] 48.74725 51.25613 47.96276
##  [5,] 50.20091 51.78332 49.39795
##  [6,] 59.77977 56.80869 55.84425
##  [7,] 53.47260 52.57011 57.50065
##  [8,] 70.80699 68.55886 62.79705
##  [9,] 66.69002 68.30327 65.35426
## [10,] 60.65127 61.74731 68.94484
# plotting the following price-paths: minimum ST, 25%, 50%, 75% and maximum St
plot(st[1,], type="l", xlab="Steps", ylab="Prices", 
     ylim=c(st[1, n+1]*.98,st[m, n+1]*1.02), col="red")
lines(st[round(m*.25),], col="orange")
lines(st[round(m*.50),], col="grey")
lines(st[round(m*.75),], col="blue")
lines(st[m,], col="green")

Custom Function

We can define a custom function so that we can apply the method to a variety of options quickly.

For ease of use, I include default values for all inputs in this didactic example, but in a real setting I would only include a default value for the number of simulations m and the number of steps n.

It wouldn’t make sense to include default values for the other inputs since in those cases I prefer to get an error message than to make a wrong assumption.

Additionally, I’ll leave the price-paths graph and the table with the prices for each option type, for didactic purposes, but this function should only calculate the required option type.

mca <- function(m = 100, n = 10, s = 52, k = 50, 
                t = 0.50, v = 0.20, r = 0.10) {
  dt <- t/n
  st <- matrix(s, m, n+1)
  for (j in 2:(n+1)) {
    epsilon <- rnorm(m)
    st[, j] <- st[, j-1] * exp((r - (v^2)/2)*dt + v*epsilon*sqrt(dt))  
  }
  stavg <- rowMeans(st)
  
  st <- st[order(st[, n+1]),]
  
  plot(st[1,], type="l", xlab="Steps", ylab="Prices",
       ylim=c(st[1,n+1]*.98,st[m,n+1]*1.02), col="red")
  lines(st[round(m*.25),], col="orange")
  lines(st[round(m*.50),], col="grey")
  lines(st[round(m*.75),], col="blue")
  lines(st[m,], col="green")

  op <- matrix(NA,2,3)
  rownames(op) <- c("call","put")
  colnames(op) <- c("eur","afix","afloat")
  op[1, 1] <- mean(pmax(st[, n+1] - k, 0))
  op[2, 1] <- mean(pmax(k - st[, n+1], 0))
  op[1, 2] <- mean(pmax(stavg - k, 0))
  op[2, 2] <- mean(pmax(k - stavg, 0))
  op[1, 3] <- mean(pmax(st[, n+1] - stavg, 0))
  op[2, 3] <- mean(pmax(stavg - st[, n+1], 0))
  
  op <- op*exp(-r*t)
  
  # For a large m and n, we don't want to display the whole tree,
  # so we want only the final option value, where is it stored?
  return(___) 
}
# To verify it was created successfully
mca
## function(m = 100, n = 10, s = 52, k = 50, 
##                 t = 0.50, v = 0.20, r = 0.10) {
##   dt <- t/n
##   st <- matrix(s, m, n+1)
##   for (j in 2:(n+1)) {
##     epsilon <- rnorm(m)
##     st[, j] <- st[, j-1] * exp((r - (v^2)/2)*dt + v*epsilon*sqrt(dt))  
##   }
##   stavg <- rowMeans(st)
##   
##   st <- st[order(st[, n+1]),]
##   
##   plot(st[1,], type="l", xlab="Steps", ylab="Prices",
##        ylim=c(st[1,n+1]*.98,st[m,n+1]*1.02), col="red")
##   lines(st[round(m*.25),], col="orange")
##   lines(st[round(m*.50),], col="grey")
##   lines(st[round(m*.75),], col="blue")
##   lines(st[m,], col="green")
## 
##   op <- matrix(NA,2,3)
##   rownames(op) <- c("call","put")
##   colnames(op) <- c("eur","afix","afloat")
##   op[1, 1] <- mean(pmax(st[, n+1] - k, 0))
##   op[2, 1] <- mean(pmax(k - st[, n+1], 0))
##   op[1, 2] <- mean(pmax(stavg - k, 0))
##   op[2, 2] <- mean(pmax(k - stavg, 0))
##   op[1, 3] <- mean(pmax(st[, n+1] - stavg, 0))
##   op[2, 3] <- mean(pmax(stavg - st[, n+1], 0))
##   
##   op <- op*exp(-r*t)
##   
##   # For a large m and n, we don't want to display the whole tree,
##   # so we want only the final option value, where is it stored?
##   return(op[1,1])
## }
# Function use with default values
mca()

## [1] 5.909553
# Test the function with your own values, for example:
mca(m = 10000, n = 1000)

## [1] 5.482699