Homework 8- European Option Pricing

Instructor: Le Nhat Tan

Libraries

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(ggplot2)
library(VNDS) # https://github.com/phamdinhkhanh/VNDS
library(yuima)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: stats4

## Loading required package: expm

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 
## Attaching package: 'expm'

## The following object is masked from 'package:Matrix':
## 
##     expm

## Loading required package: cubature

## Loading required package: mvtnorm

## ########################################

## This is YUIMA Project package v.1.15.15

## Why don't you try yuimaGUI package?

## Visit: http://www.yuima-project.com

## ########################################

## 
## Attaching package: 'yuima'

## The following object is masked from 'package:tibble':
## 
##     char

## The following object is masked from 'package:stats':
## 
##     simulate

Inclass

HAG

HAG = tq_get(symbol = 'HAG', from = '2018-01-01', to = '2022-09-30',
             src = 'CAFEF')

## #HAG from 2018-01-01 to 2022-09-30 already cloned

HAG = as.data.frame(HAG)
plot(HAG$date, HAG$adjusted, type = 'l', col = 'blue',
     xlab = 'Time', ylab = 'Adjusted Price', main = 'HAG')

price = rev(HAG$adjusted)
delta = 1 / 252
gBm = setModel(drift = "mu*x", diffusion = "sigma*x")

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

mod = setYuima(model = gBm, data = setData(price, delta = delta))
fit = qmle(mod, start = list(mu = 1, sigma = 1),
           lower = list(mu = 0.001, sigma = 0.001),
           upper = list(mu = 10, sigma = 10))
coef(summary(fit))

##        Estimate  Std. Error
## sigma 0.4770663 0.009879862
## mu    0.2407408 0.220184436

Estimated parameters are \(\mu=0.2407\) and \(\sigma=0.4771\)

tail(HAG, 1)

##            date open high  low close  volume adjusted
## 1184 2018-01-02 7.35 7.38 7.23  7.25 5924880     7.25

(length(price) - 1) / 252

## [1] 4.694444

geometric_brownian_motion = function(start, time, size, mu, 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 = mu * dt + sigma * dWt
  return(cumprod(c(start, dSt_by_St + 1)))
}

sigma = coef(summary(fit))[1]; mu = coef(summary(fit))[2]
start = price[length(price)]; x = HAG$date[1] + 1:365
time = 1; size = length(x) - 1; num_path = 500; val = NA;
for (i in 1:num_path) {
  y = geometric_brownian_motion(start, time, size, mu, sigma)
  val = rbind(val, data.frame(x, y, path = i, col = 'black'))
}

val_HAG = val %>% group_by(x) %>% summarize(y = mean(y))
val_HAG = rbind(val_HAG, data.frame(x = HAG$date, y = HAG$adjusted))
val_HAG$path = num_path + 1
val_HAG$col = 'red'
val = rbind(val, val_HAG)

ggplot(val, aes(x, y, group = path)) +
  geom_line(data = val[val$col == 'black',], aes(color = 'GBM Sample Paths'),
            alpha = 0.2) +
  geom_line(data = val[val$col == 'red',], aes(color = 'HAG as Average'), alpha = 2,
            size = 1.1) +
  scale_color_manual(values = c('black', 'red')) +
  theme_classic() +
  theme(legend.title = element_blank()) +
  xlab('Time') +
  ylab('Price') +
  ggtitle('Forecasting HAG Future movements')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 2 row(s) containing missing values (geom_path).

## VIC

VIC = tq_get(symbol = 'VIC', from = '2018-01-01', to = '2022-09-30',
             src = 'CAFEF')

## #VIC from 2018-01-01 to 2022-09-30 already cloned

VIC = as.data.frame(VIC)
plot(VIC$date, VIC$adjusted, type = 'l', col = 'blue',
     xlab = 'Time', ylab = 'Adjusted Price', main = 'VIC')

price = rev(VIC$adjusted)
gBm = setModel(drift = "mu*x", diffusion = "sigma*x")

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

mod = setYuima(model = gBm, data = setData(price, delta = delta))
fit = qmle(mod, start = list(mu = 1, sigma = 1),
           lower = list(mu = 0.001, sigma = 0.001),
           upper = list(mu = 10, sigma = 10))
coef(summary(fit))

##         Estimate Std. Error
## sigma 0.30046910 0.00637214
## mu    0.03401463 0.13867805

Estimated parameters are \(\mu=0.03401463\) and \(\sigma=0.30046910\)

tail(VIC, 1)

##            date open high  low close  volume adjusted
## 1184 2018-01-02 77.6 78.5 77.4  78.2 1435130    57.45

(length(price) - 1) / 252

## [1] 4.694444

sigma = coef(summary(fit))[1]; mu = coef(summary(fit))[2]
start = price[length(price)]; x = VIC$date[1] + 1:365
size = length(x) - 1; val = NA;
for (i in 1:num_path) {
  y = geometric_brownian_motion(start, time, size, mu, sigma)
  val = rbind(val, data.frame(x, y, path = i, col = 'black'))
}

val_VIC = val %>% group_by(x) %>% summarize(y = mean(y))
val_VIC = rbind(val_VIC, data.frame(x = VIC$date, y = VIC$adjusted))
val_VIC$path = num_path + 1
val_VIC$col = 'red'
val = rbind(val, val_VIC)

ggplot(val, aes(x, y, group = path)) +
  geom_line(data = val[val$col == 'black',],
            aes(color = 'GBM Sample Paths'), alpha = 0.2) +
  geom_line(data = val[val$col == 'red',], aes(color = 'VIC'),
            alpha = 2, size = 1.1) +
  scale_color_manual(values = c('black', 'red')) +
  theme_classic() +
  theme(legend.title = element_blank()) +
  xlab('Time') +
  ylab('Price') +
  ggtitle('Forecasting VIC Future movements')

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 2 row(s) containing missing values (geom_path).

Homework

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) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  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) {
  tree = matrix(0, nrow = nrow(stock_tree), ncol = nrow(stock_tree))
  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)
}

Problem 2

sigma = 0.3; k = 8; r = 0.07
time = 1 - 0.5; price = 10
step = 500
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(price, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_call(st, p, r, delta_t, k)[1][1]

## [1] 2.375834

Problem 3

k = 10; time = 1 - 0.5
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(price, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, k)[1][1]

## [1] 0.6690108

Problem 5

We assume the time period January 1 2022 − September 30 2022.

stock_list = c('HPG', 'VHM', 'SSI', 'MBB', 'MWG')
for (stock in stock_list) {
    data = tq_get(symbol = stock, from = '2022-01-01', to = '2022-09-30',
                  src = 'CAFEF')
    data = as.data.frame(data)
    price = rev(data$adjusted)
    delta = 1 / 252
    gBm = setModel(drift = "mu*x", diffusion = "sigma*x")
    mod = setYuima(model = gBm, data = setData(price, delta = delta))
    fit = qmle(mod, start = list(mu = 1, sigma = 1),
               lower = list(mu = 0.001, sigma = 0.001),
               upper = list(mu = 10, sigma = 10))
    sigma = coef(summary(fit))[1]; mu = coef(summary(fit))[2]
    print(paste0('Stock: ', stock, ', mu: ', mu, ', sigma: ', sigma))
    print('')
}

## #HPG from 2022-01-01 to 2022-09-30 already cloned

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

## [1] "Stock: HPG, mu: 0.00105574899140394, sigma: 0.38493408076749"
## [1] ""
## #VHM from 2022-01-01 to 2022-09-30 already cloned

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

## [1] "Stock: VHM, mu: 0.00105574899140394, sigma: 0.275376909394133"
## [1] ""
## #SSI from 2022-01-01 to 2022-09-30 already cloned

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

## [1] "Stock: SSI, mu: 0.00105574899140394, sigma: 0.512427009892711"
## [1] ""
## #MBB from 2022-01-01 to 2022-09-30 already cloned

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

## [1] "Stock: MBB, mu: 0.00105574899140394, sigma: 0.357350054322151"
## [1] ""
## #MWG from 2022-01-01 to 2022-09-30 already cloned

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

## [1] "Stock: MWG, mu: 0.00105574899140394, sigma: 0.365469521928479"
## [1] ""

Problem 6

Binomial Pricing

VCB = tq_get(symbol = 'VCB', from = '2014-02-01', to = '2019-02-01',
             src = 'CAFEF')

## #VCB from 2014-02-01 to 2019-02-01 already cloned

sigma = sd(VCB$close) / sqrt(5)
sigma

## [1] 5.146645

k = 62; r = 0.08; time = 0.5
VCB = as.data.frame(VCB)
price = rev(VCB$close)
price = price[length(price)]
delta_t = time / step
u = exp(sigma * sqrt(delta_t))
d = exp(-sigma * sqrt(delta_t))
st = stock_tree(price, u, d, step + 1)
p = (exp(r * delta_t) - d) / (u - d)
European_put(st, p, r, delta_t, k)[1][1]

## [1] 55.57305

Black-Scholes-Merton

price = rev(VCB$close)
delta = 1 / 252
gBm = setModel(drift = "mu*x", diffusion = "sigma*x")

## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.

mod = setYuima(model = gBm, data = setData(price, delta = delta))
fit = qmle(mod, start = list(mu = 1, sigma = 1),
           lower = list(mu = 0.001, sigma = 0.001),
           upper = list(mu = 10, sigma = 10))
sigma = coef(summary(fit))[1]; mu = coef(summary(fit))[2]
print(paste0('mu: ', mu, ', sigma: ', sigma))

## [1] "mu: 0.200453262009897, sigma: 0.340860249618946"

sigma = sigma * sqrt(252)
price = price[length(price)]
d_plus = (log(price) - log(k) +
            (r + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d_minus = (log(price) - log(k) +
             (r - sigma * sigma / 2) * time) / (sigma * sqrt(time))
k * exp(-r * time) * pnorm(-d_minus) - price * pnorm(-d_plus)

## [1] 56.33554