Financial Mathematics 2
Numerical Methods

Instructor: Dr. Le Nhat Tan

1 Libraries

library(tidyverse)
library(plyr)

2 Binomial Tree

Note that the binomial tree method only applies to the geometric Brownian motion \[dS_t=rS_tdt+\sigma S_tdW_t^{\mathbb{Q}}\]

stock_tree = function(stock, up, down, n) {
  tree = matrix(0, nrow = n, ncol = n)
  for (i in 1:n) {
    for (j in 1:i) {
      tree[i,j] = stock * up^(j - 1) * down^((i - 1) - (j - 1))
    }
  }
  return(tree)
}

European_call = function(stock_tree, p, r, delta_t, strike, digital = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (digital) {
    tree[nrow(tree),] = (stock_tree[nrow(stock_tree),] > strike) * 1
  } else {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      tree[i, j] = ((1 - p) * tree[i + 1, j] + 
                      p * tree[i + 1, j + 1])/exp(r * delta_t)
    }
  }
  return(tree)
}

European_put = function(stock_tree, p, r, delta_t, strike, digital = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (digital) {
    tree[nrow(tree),] = (stock_tree[nrow(stock_tree),] < strike) * 1
  } else {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      tree[i, j] = ((1 - p) * tree[i + 1, j] + 
                      p * tree[i + 1, j + 1])/exp(r * delta_t)
    }
  }
  return(tree)
}

American_call = function(stock_tree, p, r, delta_t, strike, digital = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (digital) {
    tree[nrow(tree),] = (stock_tree[nrow(stock_tree),] > strike) * 1
  } else {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      no_exercise = ((1 - p) * tree[i + 1, j] +
                       p * tree[i + 1, j + 1])/exp(r * delta_t)
      
      if (digital) {
        tree[i, j] = (stock_tree[i, j] >= strike) * 1 +
          (stock_tree[i, j] < strike) * no_exercise
      } else {
        tree[i, j] = max(stock_tree[i, j] - strike, no_exercise)
      }
    }
  }
  return(tree)
}

American_put = function(stock_tree, p, r, delta_t, strike, digital = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (digital) {
    tree[nrow(tree),] = (stock_tree[nrow(stock_tree),] < strike) * 1
  } else {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      no_exercise = ((1 - p) * tree[i + 1, j] +
                       p * tree[i + 1, j + 1])/exp(r * delta_t)
      if (digital) {
        tree[i, j] = (stock_tree[i, j] < strike) * 1 +
          (stock_tree[i, j] >= strike) * no_exercise
      } else {
        tree[i, j] = max(strike - stock_tree[i, j], no_exercise)
      }
    }
  }
  return(tree)
}

European_call_out = function(stock_tree, p, r, delta_t, strike, barrier, up = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (up) {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0) * (stock_tree[nrow(stock_tree),] < barrier)
  } else {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0) * (stock_tree[nrow(stock_tree),] > barrier)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      if ((up & stock_tree[i, j] > barrier) | ((!up) & stock_tree[i, j] < barrier)) {tree[i, j] = 0}
      else {
        tree[i, j] = ((1 - p) * tree[i + 1, j] + 
                      p * tree[i + 1, j + 1])/exp(r * delta_t)
      }
    }
  }
  return(tree)
}

European_put_out = function(stock_tree, p, r, delta_t, strike, barrier, up = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (up) {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0) * (stock_tree[nrow(stock_tree),] < barrier)
  } else {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0) * (stock_tree[nrow(stock_tree),] > barrier)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      if ((up & stock_tree[i, j] > barrier) | ((!up) & stock_tree[i, j] < barrier)) {tree[i, j] = 0}
      else {
        tree[i, j] = ((1 - p) * tree[i + 1, j] + 
                      p * tree[i + 1, j + 1])/exp(r * delta_t)
      }
    }
  }
  return(tree)
}

American_call_out = function(stock_tree, p, r, delta_t, strike, barrier, up = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (up) {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0) * (stock_tree[nrow(stock_tree),] < barrier)
  } else {
    tree[nrow(tree),] = pmax(stock_tree[nrow(stock_tree),] - strike, 0) * (stock_tree[nrow(stock_tree),] > barrier)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      if ((up & stock_tree[i, j] > barrier) | ((!up) & stock_tree[i, j] < barrier)) {tree[i, j] = 0}
      else {
        no_exercise = ((1 - p) * tree[i + 1, j] +
                       p * tree[i + 1, j + 1])/exp(r * delta_t)
        tree[i, j] = max(stock_tree[i, j] - strike, no_exercise)
      }
    }
  }
  return(tree)
}

American_put_out = function(stock_tree, p, r, delta_t, strike, barrier, up = F) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  if (up) {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0) * (stock_tree[nrow(stock_tree),] < barrier)
  } else {
    tree[nrow(tree),] = pmax(strike - stock_tree[nrow(stock_tree),], 0) * (stock_tree[nrow(stock_tree),] > barrier)
  }
  for (i in (nrow(tree) - 1):1) {
    for (j in 1:i) {
      if ((up & stock_tree[i, j] > barrier) | ((!up) & stock_tree[i, j] < barrier)) {tree[i, j] = 0}
      else {
        no_exercise = ((1 - p) * tree[i + 1, j] +
                       p * tree[i + 1, j + 1])/exp(r * delta_t)
        tree[i, j] = max(strike - stock_tree[i, j], no_exercise)
      }
    }
  }
  return(tree)
}

2.1 Call Options

sigma = 0.3; K = 10; r = 0.1; t = 0.5; T = 1; St = 9
time = T - t
# Vanilla Call - Analytical (Black-Scholes)
d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
St * pnorm(d_plus) - exp(-r * time) * K * pnorm(d_minus)

## [1] 0.5520873

# Vanilla Call - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
st

##          [,1]      [,2]     [,3]     [,4]     [,5]
## [1,] 9.000000  0.000000  0.00000  0.00000  0.00000
## [2,] 8.094287 10.007058  0.00000  0.00000  0.00000
## [3,] 7.279721  9.000000 11.12680  0.00000  0.00000
## [4,] 6.547128  8.094287 10.00706 12.37184  0.00000
## [5,] 5.888260  7.279721  9.00000 11.12680 13.75619

p = (exp(r * delta_t) - d) / (u - d)
European_call(st, p, r, delta_t, K)

##           [,1]      [,2]      [,3]     [,4]     [,5]
## [1,] 0.5905478 0.0000000 0.0000000 0.000000 0.000000
## [2,] 0.1640556 0.9786348 0.0000000 0.000000 0.000000
## [3,] 0.0000000 0.3118480 1.5866853 0.000000 0.000000
## [4,] 0.0000000 0.0000000 0.5927818 2.496058 0.000000
## [5,] 0.0000000 0.0000000 0.0000000 1.126800 3.756186

# Vanilla Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call(st, p, r, delta_t, K)[1, 1]

## [1] 0.5520472

# American Call - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call(st, p, r, delta_t, K)

##           [,1]      [,2]      [,3]     [,4]     [,5]
## [1,] 0.5905478 0.0000000 0.0000000 0.000000 0.000000
## [2,] 0.1640556 0.9786348 0.0000000 0.000000 0.000000
## [3,] 0.0000000 0.3118480 1.5866853 0.000000 0.000000
## [4,] 0.0000000 0.0000000 0.5927818 2.496058 0.000000
## [5,] 0.0000000 0.0000000 0.0000000 1.126800 3.756186

# American Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call(st, p, r, delta_t, K)[1, 1]

## [1] 0.5520472

# Digital Call - Analytical
exp(-r * time) * pnorm(d_minus)

## [1] 0.3393942

# Digital Call - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call(st, p, r, delta_t, K, digital = T)

##           [,1]      [,2]      [,3]      [,4] [,5]
## [1,] 0.3453618 0.0000000 0.0000000 0.0000000    0
## [2,] 0.1455942 0.5287640 0.0000000 0.0000000    0
## [3,] 0.0000000 0.2767554 0.7623255 0.0000000    0
## [4,] 0.0000000 0.0000000 0.5260755 0.9875778    0
## [5,] 0.0000000 0.0000000 0.0000000 1.0000000    1

# Digital Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call(st, p, r, delta_t, K, digital = T)[1, 1]

## [1] 0.3473722

# American Digital Call - Analytical
lambda_plus = (-(r - sigma ^ 2 / 2) +
                 sqrt((r - sigma ^ 2 / 2) ^ 2 +
                        2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = (-(r - sigma ^ 2 / 2) -
                  sqrt((r - sigma ^ 2 / 2) ^ 2 +
                         2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (log(St) - log(K) +
            time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                          2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) -
             time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                           2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
(St / K) ^ lambda_plus * pnorm(d_plus) +
  (St / K) ^ lambda_minus * pnorm(d_minus)

## [1] 0.648492

# American Digital Call - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call(st, p, r, delta_t, K, digital = T)

##           [,1]      [,2] [,3] [,4] [,5]
## [1,] 0.6537987 0.0000000    0    0    0
## [2,] 0.2767554 1.0000000    0    0    0
## [3,] 0.0000000 0.5260755    1    0    0
## [4,] 0.0000000 0.0000000    1    1    0
## [5,] 0.0000000 0.0000000    0    1    1

# American Digital Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call(st, p, r, delta_t, K, digital = T)[1, 1]

## [1] 0.6424439

2.2 Put Options

sigma = 0.3; K = 25; r = 0.07; t = 0; T = 1.5; St = 30
time = T - t
# Vanilla Put - Analytical (Black-Scholes)
d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
K * exp(-r * time) * pnorm(-d_minus) - St * pnorm(-d_plus)

## [1] 1.172938

# Vanilla Put - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
st

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 30.00000  0.00000  0.00000  0.00000  0.00000
## [2,] 24.96527 36.05008  0.00000  0.00000  0.00000
## [3,] 20.77549 30.00000 43.32028  0.00000  0.00000
## [4,] 17.28886 24.96527 36.05008 52.05665  0.00000
## [5,] 14.38737 20.77549 30.00000 43.32028 62.55488

p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, K)

##           [,1]      [,2] [,3] [,4] [,5]
## [1,]  1.333102 0.0000000    0    0    0
## [2,]  2.427134 0.4153394    0    0    0
## [3,]  4.259416 0.8998999    0    0    0
## [4,]  7.063431 1.9497785    0    0    0
## [5,] 10.612634 4.2245101    0    0    0

# Vanilla Put - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, K)[1, 1]

## [1] 1.172921

# American Put - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_put(st, p, r, delta_t, K)

##           [,1]      [,2] [,3] [,4] [,5]
## [1,]  1.396782 0.0000000    0    0    0
## [2,]  2.565109 0.4153394    0    0    0
## [3,]  4.558360 0.8998999    0    0    0
## [4,]  7.711143 1.9497785    0    0    0
## [5,] 10.612634 4.2245101    0    0    0

# American Put - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_put(st, p, r, delta_t, K)[1, 1]

## [1] 1.286398

# Digital Put - Analytical
exp(-r * time) * pnorm(-d_minus)

## [1] 0.2474335

# Digital Put - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, K, digital = T)

##           [,1]       [,2] [,3] [,4] [,5]
## [1,] 0.2469464 0.00000000    0    0    0
## [2,] 0.4258661 0.09831658    0    0    0
## [3,] 0.6861448 0.21301876    0    0    0
## [4,] 0.9740915 0.46153955    0    0    0
## [5,] 1.0000000 1.00000000    0    0    0

# Digital Put - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, K, digital = T)[1, 1]

## [1] 0.2540443

# American Digital Put - Analytical
lambda_plus = ((r - sigma ^ 2 / 2) +
                 sqrt((r - sigma ^ 2 / 2) ^ 2 +
                        2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = ((r - sigma ^ 2 / 2) -
                  sqrt((r - sigma ^ 2 / 2) ^ 2 +
                         2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (log(K) - log(St) +
            time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                          2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
d_minus = (log(K) - log(St) -
             time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                           2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
(K / St) ^ lambda_plus * pnorm(d_plus) +
  (K / St) ^ lambda_minus * pnorm(d_minus)

## [1] 0.5690159

# American Digital Put - Binomial Tree (4 Step)
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_put(st, p, r, delta_t, K, digital = T)

##           [,1]      [,2] [,3] [,4] [,5]
## [1,] 0.5707227 0.0000000    0    0    0
## [2,] 1.0000000 0.2130188    0    0    0
## [3,] 1.0000000 0.4615396    0    0    0
## [4,] 1.0000000 1.0000000    0    0    0
## [5,] 1.0000000 1.0000000    0    0    0

# American Digital Put - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_put(st, p, r, delta_t, K, digital = T)[1, 1]

## [1] 0.5620131

2.3 Barrier Options

# European Down-and-Out Call - Analytical
sigma = 0.2206; K = 45; r = 0.07; t = 0; T = 0.5; St = 42.75; barrier = 38
time = T - t
vanilla_call = function(St) {
  d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
  d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
  return(St * pnorm(d_plus) - exp(-r * time) * K * pnorm(d_minus))
}

lambda = 1 - 2 * r / sigma^2
vanilla_call(St) - (St / barrier)^lambda * vanilla_call(barrier^2 / St)

## [1] 2.235596

# European Down-and-Out Call - Binomial Tree (4 Step)

step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call_out(st, p, r, delta_t, K, barrier)

##           [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 2.4489052 0.000000 0.000000 0.000000  0.00000
## [2,] 0.7483295 3.956442 0.000000 0.000000  0.00000
## [3,] 0.0000000 1.406322 6.221741 0.000000  0.00000
## [4,] 0.0000000 0.000000 2.642876 9.411842  0.00000
## [5,] 0.0000000 0.000000 0.000000 4.966708 13.40168

# European Down-and-Out Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call_out(st, p, r, delta_t, K, barrier)[1, 1]

## [1] 2.237978

# American Down-and-Out Call - Binomial Tree (4 Step)
sigma = 0.2; K = 100; r = 0.05; t = 0; T = 1; St = 90; barrier = 80
time = T - t
step = 4
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call_out(st, p, r, delta_t, K, barrier)

##          [,1]     [,2]      [,3]      [,4]     [,5]
## [1,] 5.442102 0.000000  0.000000  0.000000  0.00000
## [2,] 1.487242 8.968184  0.000000  0.000000  0.00000
## [3,] 0.000000 2.800160 14.478725  0.000000  0.00000
## [4,] 0.000000 0.000000  5.272104 22.729513  0.00000
## [5,] 0.000000 0.000000  0.000000  9.926248 34.26422

# American Down-and-Out Call - Binomial Tree (1000 Step)
step = 1000
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(St, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
American_call_out(st, p, r, delta_t, K, barrier)[1, 1]

## [1] 4.672202

3 Monte Carlo

3.1 GBM

\[\begin{align*} S_T &= S_t\cdot\exp\left(\left(r-\frac{\sigma^2}{2}\right)(T-t)+\sigma\left(W_T^{\mathbb{Q}}-W_t^{\mathbb{Q}}\right)\right)\\ &\sim S_t\cdot\exp\left(\mathcal{N}\left(\left(r-\frac{\sigma^2}{2}\right)(T-t), \sigma^2(T-t)\right)\right) \end{align*}\]

3.1.1 Call Options

geometric_brownian_motion = function(S0, time, size, r, sigma) {
  # dSt = St(\mu dt + \sigma dWt)
  step = time / size
  dWt = rnorm(size, mean = 0, sd = sqrt(step))
  dt = rep(step, length(dWt))
  dSt_by_St = r * dt + sigma * dWt
  return(cumprod(c(S0, dSt_by_St + 1)))
}

american_digital = function(S, r, K, time, call = T) {
  dt = time / (length(S) - 1)
  price = 0
  for (j in 1:length(S)){
    if ((S[j] >= K & call) | (S[j] <= K & !call)){
      price = exp(-r * (j - 1) * dt)
      break
    }
  }
  return(price)
}

american_digital_GBM = function(St, t, T, r, sigma, K,
                                call = T, size = 1000, rep = 1000) {
  time = T - t
  stocks = replicate(rep, geometric_brownian_motion(St, time, size, r, sigma))
  df = data.frame(stocks) %>%
    (colwise(function(x) {return(american_digital(x, r, K, time, call))}))
  return(mean(as.matrix(df)))
}

european_out = function(S, r, K, time, barrier, up = F, call = T) {
  if ((up & (max(S) > barrier)) | ((!up) & (min(S) < barrier))) {return(0)}
  ST = S[length(S)]
  if (call) {return((ST > K) * (ST - K) * exp(-r * time))}
  return((ST < K) * (K - ST) * exp(-r * time))
}

european_out_GBM = function(St, t, T, r, sigma, K, barrier,
                            up = F, call = T, size = 1000, rep = 1000) {
  time = T - t
  stocks = replicate(rep, geometric_brownian_motion(St, time, size, r, sigma))
  df = data.frame(stocks) %>%
    (colwise(function(x) {return(european_out(x, r, K, time, barrier, up, call))}))
  return(mean(as.matrix(df)))
}

sigma = 0.3; K = 10; r = 0.1; t = 0.5; T = 1; St = 9
time = T - t
# Vanilla Call - Analytical (Black-Scholes)
d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
St * pnorm(d_plus) - exp(-r * time) * K * pnorm(d_minus)

## [1] 0.5520873

# Vanilla Call - Monte Carlo
sample = 100000
ST = St * exp(rnorm(sample, mean = (r - sigma ^ 2 / 2) * time,
                    sd = sigma * sqrt(time)))
mean((ST > K) * (ST - K)) * exp(-r * time)

## [1] 0.5547082

# Digital Call - Analytical
exp(-r * time) * pnorm(d_minus)

## [1] 0.3393942

# Digital Call - Monte Carlo
mean((ST > K)) * exp(-r * time)

## [1] 0.3393416

# American Digital Call - Analytical
lambda_plus = (-(r - sigma ^ 2 / 2) +
                 sqrt((r - sigma ^ 2 / 2) ^ 2 +
                        2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = (-(r - sigma ^ 2 / 2) -
                  sqrt((r - sigma ^ 2 / 2) ^ 2 +
                         2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (log(St) - log(K) +
            time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                          2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) -
             time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                           2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
(St / K) ^ lambda_plus * pnorm(d_plus) +
  (St / K) ^ lambda_minus * pnorm(d_minus)

## [1] 0.648492

# American Digital Call - Monte Carlo
american_digital_GBM(St, t, T, r, sigma, K)

## [1] 0.6402451

3.1.2 Put Options

sigma = 0.3; K = 9; r = 0.1; t = 0.5; T = 1; St = 10
time = T - t
# Vanilla Put - Analytical (Black-Scholes)
d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
K * exp(-r * time) * pnorm(-d_minus) - St * pnorm(-d_plus)

## [1] 0.2645281

# Vanilla Put - Monte Carlo
ST = St * exp(rnorm(sample, mean = (r - sigma ^ 2 / 2) * time,
                    sd = sigma * sqrt(time)))
mean((ST < K) * (K - ST)) * exp(-r * time)

## [1] 0.2629521

# Digital Put - Analytical
exp(-r * time) * pnorm(-d_minus)

## [1] 0.2526044

# Digital Put - Monte Carlo
mean((ST < K)) * exp(-r * time)

## [1] 0.2508392

# American Digital Put - Analytical
lambda_plus = ((r - sigma ^ 2 / 2) +
                 sqrt((r - sigma ^ 2 / 2) ^ 2 +
                        2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = ((r - sigma ^ 2 / 2) -
                  sqrt((r - sigma ^ 2 / 2) ^ 2 +
                         2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (log(K) - log(St) +
            time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                          2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
d_minus = (log(K) - log(St) -
             time * sqrt((r - sigma ^ 2 / 2) ^ 2 +
                           2 * sigma ^ 2 * r)) / (sigma * sqrt(time))
(K / St) ^ lambda_plus * pnorm(d_plus) +
  (K / St) ^ lambda_minus * pnorm(d_minus)

## [1] 0.5701365

# American Digital Put - Monte Carlo
american_digital_GBM(St, t, T, r, sigma, K, call = F)

## [1] 0.5766318

3.1.3 Barrier Options

# European Down-and-Out Call - Analytical
sigma = 0.2206; K = 45; r = 0.07; t = 0; T = 0.5; St = 42.75; barrier = 38
time = T - t
vanilla_call = function(St) {
  d_plus = (log(St) - log(K) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
  d_minus = (log(St) - log(K) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
  return(St * pnorm(d_plus) - exp(-r * time) * K * pnorm(d_minus))
}

lambda = 1 - 2 * r / sigma^2
vanilla_call(St) - (St / barrier)^lambda * vanilla_call(barrier^2 / St)

## [1] 2.235596

# European Down-and-Out Call - Monte Carlo
european_out_GBM(St, t, T, r, sigma, K, barrier)

## [1] 2.217315

3.2 Risk-Neutral ABM

\[\begin{align*} S_T &= S_t+r(T-t)+\sigma\left(W_T^{\mathbb{Q}}-W_t^{\mathbb{Q}}\right)\\ &\sim \mathcal{N}(S_t+r(T-t),\sigma^2(T-t)) \end{align*}\]

3.2.1 Call Options

drifted_brownian_motion = function(start, time, size, mu, sigma) {
  # dXt = mu dt + sigma dWt
  step = time / size
  dWt = rnorm(size, mean = 0, sd = sqrt(step))
  dt = rep(step, length(dWt))
  dXt = mu * dt + sigma * dWt
  return(cumsum(c(start, dXt)))
}

american_digital_ABM = function(St, t, T, r, sigma, K,
                                call = T, size = 1000, rep = 1000) {
  time = T - t
  stocks = replicate(rep, drifted_brownian_motion(St, time, size, r, sigma))
  df = data.frame(stocks) %>%
    (colwise(function(x) {
        if (call) {
          return(american_digital(x, r, K, time))
        }
        return(american_digital(x, r, K, time, call = F))
      }))
  return(mean(as.matrix(df)))
}

sigma = 0.3; K = 10; r = 0.1; t = 0.5; T = 1; St = 9.9
time = T - t
# Vanilla Call - Analytical (Black-Scholes)
sigma_hat = exp(-r * time) * sigma * sqrt(time)
d = (St + r * time - K) / (sigma * sqrt(time))
sigma_hat * d * pnorm(d) + sigma_hat * dnorm(d)

## [1] 0.05894617

# Vanilla Call - Monte Carlo
sample = 100000
ST = rnorm(sample, mean = St + r * time,
                    sd = sigma * sqrt(time))
mean((ST > K) * (ST - K)) * exp(-r * time)

## [1] 0.05908891

# Digital Call - Analytical
exp(-r * time) * pnorm(d)

## [1] 0.3869904

# Digital Call - Monte Carlo
mean((ST > K)) * exp(-r * time)

## [1] 0.3892621

# American Digital Call - Analytical
lambda_plus = (-r + sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = (-r - sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (St - K + (T - t) *
            sqrt(r^2 + 2 * sigma ^ 2 * r)) / (sigma * sqrt(T - t))
d_minus = (St - K - (T - t) *
             sqrt(r^2 + 2 * sigma ^ 2 * r)) / (sigma * sqrt(T - t))
exp(lambda_plus * (St - K)) * pnorm(d_plus) +
  exp(lambda_minus * (St - K)) * pnorm(d_minus)

## [1] 0.6956341

# American Digital Call - Monte Carlo
american_digital_ABM(St, t, T, r, sigma, K)

## [1] 0.6597976

3.2.2 Put Options

sigma = 0.3; K = 9.9; r = 0.01; t = 0.5; T = 1; St = 10
time = T - t
# Vanilla Put - Analytical (Black-Scholes)
sigma_hat = exp(-r * time) * sigma * sqrt(time)
d = (St + r * time - K) / (sigma * sqrt(time))
sigma_hat * dnorm(-d) - sigma_hat * d * pnorm(-d)

## [1] 0.04207792

# Vanilla Put - Monte Carlo
ST = rnorm(sample, mean = St + r * time,
                    sd = sigma * sqrt(time))
mean((ST < K) * (K - ST))

## [1] 0.04208951

# Digital Put - Analytical
exp(-r * time) * pnorm(-d)

## [1] 0.3087613

# Digital Put - Monte Carlo
mean((ST < K)) * exp(-r * time)

## [1] 0.3087922

# American Digital Put - Analytical
lambda_plus = (r + sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma ^ 2)
lambda_minus = (r - sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma ^ 2)
d_plus = (K - St + (T - t) *
            sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma * sqrt(T - t))
d_minus = (K - St - (T - t) *
             sqrt(r ^ 2 + 2 * sigma ^ 2 * r)) / (sigma * sqrt(T - t))
exp(lambda_plus * (K - St)) * pnorm(d_plus) +
  exp(lambda_minus * (K - St)) * pnorm(d_minus)

## [1] 0.6292933

# American Digital Put - Monte Carlo
american_digital_ABM(St, t, T, r, sigma, K, call = F)

## [1] 0.6139867

3.3 Real-World ABM

\[\begin{align*} S_T &= e^{r(T-t)}S_t+\sigma\int_t^Te^{r(T-s)}dW_s^{\mathbb{Q}}\\ &\sim \mathcal{N}\left(e^{r(T-t)}S_t,\frac{\sigma^2\left(e^{2r(T-t)}-1\right)}{2r}\right) \end{align*}\]

3.3.1 Call Options

sigma = 0.3; K = 10; r = 0.1; t = 0.5; T = 1; St = 9.9
time = T - t
# Vanilla Call - Analytical (Black-Scholes)
sigma_hat = sigma * sqrt(time)
d = (St + r * time - K) / (sigma * sqrt(time))
sigma_hat * d * pnorm(d) + sigma_hat * dnorm(d)

## [1] 0.0619684

# Vanilla Call - Monte Carlo
sample = 100000
ST = rnorm(sample, mean = St + r * time,
                    sd = sigma * sqrt(time))
mean((ST > K) * (ST - K)) * exp(-r * time)

## [1] 0.05922303

# Digital Call - Analytical
exp(-r * time) * pnorm(d)

## [1] 0.3869904

# Digital Call - Monte Carlo
mean((ST > K)) * exp(-r * time)

## [1] 0.3902133

3.3.2 Put Options

sigma = 0.3; K = 25; r = 0.07; t = 0; T = 1.5; St = 30
time = T - t
# Vanilla Put - Analytical (Black-Scholes)
sigma_hat = exp(-r * time) * sigma * sqrt(time)
d = (St + r * time - K) / (sigma * sqrt(time))
sigma_hat * dnorm(-d) - sigma_hat * d * pnorm(-d)

## [1] 8.111083e-46

# Vanilla Put - Monte Carlo
ST = rnorm(sample, mean = St + r * time,
                    sd = sigma * sqrt(time))
mean((ST < K) * (K - ST)) * exp(-r * time)

## [1] 0

# Digital Put - Analytical
exp(-r * time) * pnorm(-d)

## [1] 3.098491e-44

# Digital Put - Monte Carlo
mean((ST < K)) * exp(-r * time)

## [1] 0

Financial Mathematics 2 Numerical Methods