# R/LSM_american_option.R In LSMRealOptions: Value American and Real Options Through LSM Simulation

#### Documented in LSM_american_option

```#' Value American-style options through least-squares Monte Carlo (LSM) simulation
#'
#' @description
#'
#' Given a set of state variables and associated payoffs simulated through Monte Carlo simulation, solve for the value of an American-style call or put option through
#' the least-squares Monte Carlo simulation method.
#'
#' @param state_variables \code{matrix} or \code{array}. The simulated states of the underlying stochastic variables. The first dimension corresponds to the simulated values
#' of the state variables at each discrete observation point. The second dimension corresponds to each individual simulation of the state variable. The third dimension
#' corresponds to each state variable considered.
#' @param payoff \code{matrix} The payoff at each observation point resulting from exercise into the underlying asset. The first dimension corresponds to the simulated values
#' of the state variables at each discrete observation point. The second dimension corresponds to each individual simulation of the state variable.
#' @param K the exercise price of the American-style option
#' @param dt Constant, discrete time step of simulated observations
#' @param rf The annual risk-free interest rate
#' @param call \code{logical} Is the American-style option a call or put option?
#' @param orthogonal \code{character}. The orthogonal polynomial used to develop basis functions that estimate the continuation value in the LSM simulation method.
#' \code{orthogonal} arguments available are: "Power", "Laguerre", "Jacobi", "Legendre", "Chebyshev", "Hermite". See \bold{details}.
#' @param degree The number of orthogonal polynomials used in the least-squares fit. See \bold{details}.
#' @param cross_product \code{logical}. Should a cross product of state variables be considered? Relevant only when the number of state variables
#' is greater than one.
#' @param verbose \code{logical}. Should additional information be output? See \bold{values}.
#'
#' @details
#'
#'The \code{LSM_american_option} function provides an implementation of the least-squares Monte Carlo (LSM) simulation approach to numerically approximate
#'the value of American-style options (options with early exercise opportunities). The function provides flexibility in the stochastic process followed by the underlying asset, with simulated values
#'of stochastic processes provided within the \code{state_variables} argument. It also provides flexibility in the payoffs of the option, allowing for vanilla as well as more exotic options to be considered.
#'\code{LSM_american_option} also provides analysis into the exercise timing and probability of early exercise of the option.
#'
#'\bold{Least-Squares Monte Carlo Simulation:}
#'
#'The least-squares Monte Carlo (LSM) simulation method is a numeric approach first presented by Longstaff and Schwartz (2001) that
#'approximates the value of options with early exercise opportunities. The LSM simulation method is considered one of the most efficient
#'methods of valuing American-style options due to its flexibility and computational efficiency. The approach can feature multiple
#'stochastically evolving underlying uncertainties, following both standard and exotic stochastic processes.
#'
#'The LSM method first approximates stochastic variables through a stochastic process to develop cross-sectional information,
#'then directly estimates the continuation value of in-the-money simulation paths by "(regressing) the ex-post realized payoffs from
#'continuation on functions of the values of the state variables" (Longstaff and Schwartz, 2001).
#'
#'The 'LSM_american_option' function at each discrete time period, for each simulated price path, compares the payoff that results from immediate exercise of
#'the option with the expected value of continuing to hold the option for subsequent periods. The payoff of immediate exercise is provided in the \code{payoff} argument and
#'could take several different meanings depending upon the type of American-style option being valued (e.g. the current stock price, the maximum price between multiple assets, etc.).
#'
#'The immediate profit resulting from exercise of the option is dependent upon the type of option being calculated. The profit of price path \eqn{i} and time \eqn{t}
#'is given by:
#'
#'When \code{call = TRUE}:
#'\deqn{profit_{(t,i)} = max(payoff_{(t,i)} - K, 0)}{profit[t,i] = max(payoff[t,i] - K, 0)}
#'
#'When \code{call = FALSE}:
#'\deqn{profit_{(t,i)} = max(K - payoff_{(t,i)}, 0)}{profit[t,i] = max(K - payoff[t,i], 0)}
#'
#'\bold{Orthogonal Polynomials:}
#'
#'To improve the accuracy of estimation of continuation values, the economic values in each period are regressed on a linear
#'combination of a set of basis functions of the stochastic variables. These estimated regression parameters and the simulated
#'stochastic variables are then used to calculate the estimator for the expected economic values.
#'
#'Longstaff and Schwartz (2001) state that as the conditional expectation of the continuation value belongs to a Hilbert space,
#'it can be represented by a combination of orthogonal basis functions. Increasing the number of stochastic state variables
#'therefore increases the number of required basis functions exponentially. The orthogonal polynomials available in the
#'\code{LSMRealOptions} package are: Laguerre, Jacobi, Legendre (spherical), Hermite (probabilistic), Chebyshev (of the first kind).
#'The simple powers of state variables is further available. Explicit expressions of each of these orthogonal polynomials are
#'available within the textbook of Abramowitz and Stegun (1965).
#'
#' @return
#'
#'The 'LSM_american_option' function by default returns a \code{numeric} object corresponding to the calculated value of the American-style option.
#'
#'When \code{verbose = T}, 6 objects are returned within a \code{list} class object. The objects returned are:
#'
#'\tabular{ll}{
#'
#' \code{Value} \tab The calculated option value. \cr
#' \code{Standard Error} \tab The standard error of the option value. \cr
#' \code{Expected Timing} \tab The expected time of early exercise. \cr
#' \code{Expected Timing SE} \tab The standard error of the expected time of early exercise. \cr
#' \code{Exercise Probability} \tab The probability of early exercise of the option being exercised. \cr
#' \code{Cumulative Exercise Probability} \tab \code{vector}. The cumulative probability of option exercise at each discrete observation point \cr
#' }
#'
#' @references
#'
#' Abramowitz, M., and I. A. Stegun, (1965). Handbook of mathematical functions with formulas, graphs, and mathematical tables. Courier Corporation.
#'
#' Longstaff, F. A., and E. S. Schwartz, (2001). "Valuing American options by simulation: a simple least-squares approach." The review of financial
#' studies, 14(1), 113-147.
#'
#' @examples
#'
#'# Price a vanilla American put option on an asset that follows
#'# Geometric Brownian Motion
#'
#'## Step 1 - simulate stock prices:
#'stock_prices <- GBM_simulate(n = 100, t = 1, mu = 0.05,
#'                          sigma = 0.2, S0 = 100, dt = 1/2)
#'
#'## Step 2 - Value the American put option:
#'option_value <- LSM_american_option(state_variables = stock_prices,
#'                                  payoff = stock_prices,
#'                                  K = 100,
#'                                  dt = 1/2,
#'                                  rf = 0.05)
#' @export
LSM_american_option <- function(state_variables, payoff, K, dt, rf,
call = FALSE, orthogonal = "Power", degree = 2, cross_product = TRUE,
verbose = FALSE){

if(anyNA(state_variables)) stop("NA's have been specified within 'state_variables'!")

###FUNCTION DEFINITIONS:
#------------------------------------------------------------------------------
#Continuous compounding discounting
discount = function(r, t = 1) return(exp(-r*t))

combination = function(n, k){factorial(n) / (factorial(n-k)*factorial(k))}

#####orthogonal Polynomial Functions:
orthogonal_polynomials = c("Laguerre", "Jacobi", "Legendre", "Chebyshev", "Hermite")

c_m <- function(n, m, orthogonal, alpha_value, beta_value = 0){
if(orthogonal == "Laguerre")  return((-1)^m * combination(n, n-m) * (1/factorial(m)))
if(orthogonal == "Jacobi")    return(combination(n + alpha_value, m) * combination(n + beta_value, n - m))
if(orthogonal == "Legendre")  return((-1)^m * combination(n, m) * combination(2*n - 2*m, n))
if(orthogonal == "Chebyshev") return((-1)^m * factorial(n - m - 1) / (factorial(m) * factorial(n - 2*m)))
if(orthogonal == "Hermite")   return((-1)^m / (factorial(m) * (2^m) * factorial(n - 2*m)))}

g_m_x <- function(n, m, x, orthogonal){
if(orthogonal == "Laguerre") return(x^m)
if(orthogonal == "Jacobi")   return((x-1)^(n-m) * (x+1)^m)
if(orthogonal == "Legendre" || orthogonal == "Hermite") return(x^(n-2*m))
if(orthogonal == "Chebyshev") return((2*x) ^ (n-2*m))}

orthogonal_weight <- function(n, x, orthogonal, alpha_value = 0, beta_value = 0){

##If the orthogonal is "Power" then we'll just return the power:
if(orthogonal == "Power") return(x^n)

#If not, we're return an orthogonal polynomial
orthogonal_polynomials <- c("Laguerre", "Jacobi", "Legendre", "Chebyshev", "Hermite")
index <- which(orthogonal_polynomials == orthogonal)

N <- c(rep(n,2), rep(n/2, 3))[index]
if(N%%1 != 0) sprintf("warning: orthogonal weight rounded down from %i to %i", n, n-1)
N <- floor(N)
d_n <- c(1, rep(1/(2^n), 2), n/2, factorial(n))[index]

output <- 0
for(m in 0:N) output = output + c_m(n,m, orthogonal, alpha_value, beta_value) * g_m_x(n, m, x, orthogonal)

return(d_n * output)
}

###Estimated Continuation Values:
continuation_value_calc <- function(ITM, t, continuation_value, state_variables_t, orthogonal, degree, cross_product){

#We only consider the invest / delay investment decision for price paths that are in the money (ie. positive NPV):
ITM_length <- length(ITM)

#Only regress paths In the money (ITM):
if(ITM_length>0){

if(is.vector(state_variables_t)){
state_variables_t_ITM <- matrix(state_variables_t[ITM])
state_variables.ncol <- 1
} else {
state_variables.ncol <- ncol(state_variables_t)
state_variables_t_ITM <- matrix(state_variables_t[ITM,], ncol = state_variables.ncol)
}

regression.matrix <- matrix(0, nrow = ITM_length, ncol = (1 + state_variables.ncol*degree + ifelse(cross_product, factorial(state_variables.ncol - 1),0)))
index <- state_variables.ncol+1
regression.matrix[,1:index] <- c(continuation_value[ITM], state_variables_t_ITM)
index <- index + 1

##The Independent_variables includes the actual values as well as their polynomial/orthogonal weights:
Regression_data <- data.frame(matrix(0, nrow = ITM_length,
ncol = (state_variables.ncol*degree + ifelse(cross_product, factorial(state_variables.ncol - 1),0))), state_variables_t_ITM,
dependent_variable = continuation_value[ITM])

index <- 1
#Basis Functions:
for(i in 1:state_variables.ncol){
for(n in 1:degree){
Regression_data[,index] <- orthogonal_weight(n, state_variables_t_ITM[,i], orthogonal)
index <- index + 1
}
}

#Cross Product of state vectors:
if(cross_product && state_variables.ncol>1){
for(i in 1:(state_variables.ncol-1)){
for(j in (i + 1):state_variables.ncol){
Regression_data[,index] <- state_variables_t_ITM[,i] * state_variables_t_ITM[,j]
index <- index + 1
}}
}
##Finally, perform the Linear Regression
LSM_regression <- stats::lm(dependent_variable~., Regression_data)

#if(any(is.na(LSM_regression\$coefficients))) print("There are NA's present in the (complex) regression!")
##Assign fitted values
continuation_value[ITM] <- LSM_regression\$fitted.values
}
return(continuation_value)
}

if(length(K) > 1 && length(K) != G) stop("length of object 'K' does not equal 1 or number of columns of 'state_variables'!")
if(nrow(state_variables) != nrow(payoff)) stop("Dimensions of object 'state_variables' does not match 'payoff'!")

if(any(class(state_variables) == "matrix")) state_variables <- array(state_variables, dim = c(nrow(state_variables), ncol(state_variables), 1))

if(!(orthogonal %in% c("Power", "Laguerre", "Jacobi", "Legendre", "Chebyshev", "Hermite"))){
stop("Invalid orthogonal argument! 'orthogonal' must be Power, Laguerre, Jacobi, Legendre, Chebyshev or Hermite")
}

# Nominal Interest Rate:
r <- rf * dt

# Number of discrete time periods:
nperiods <- dim(state_variables)[1]
# Number of simulations:
G <- dim(state_variables)[2]

####BEGIN LSM SIMULATION

# Profit: value of immediate exercise.
profit <- matrix(0, nrow = nperiods, ncol = G)

# Continuation value: expected value of waiting to exercise. Compared with immediate profit to make exercise decisions through dynamic programming.
continuation_value <- rep(0, G)

#American option value of project, given that you can either delay or exercise:
AOV <- rep(0, G)

##Optimal period of exercise is the earliest time that exercise is triggered. If no exercise, an NA is returned:
exercise_timing <- rep(NA, G)

###Backwards induction loop
for(t in nperiods:1){
#t is representative of the index for time, but in reality the actual time period you're in is (t-1).

## STEP ONE:
### Forward foresight (high bias) - the Immediate payoff of exercising:
if(call){
profit[t,] <- pmax(payoff[t,] - K, 0)
} else {
profit[t,] <- pmax(K - payoff[t,], 0)
}

## STEP TWO:

### Estimated continuation values:
if(t < nperiods){

#We only consider the invest / delay investment decision for price paths that are in the money (ie. positive NPV):
ITM <- which(profit[t,] > 0)

#Only regress paths In the money (ITM):
if(length(ITM)>0){

continuation_value <- continuation_value_calc(ITM = ITM,
t = t,
continuation_value = AOV  * discount(r),
state_variables_t = state_variables[t,,],
orthogonal = orthogonal,
degree = degree,
cross_product = cross_product)
} else {
continuation_value <- AOV  * discount(r)
}
}

## STEP THREE:

### Dynamic programming:
exercise <- profit[t,] > continuation_value

## Discount existing values if not exercising
AOV[!exercise] <- AOV[!exercise] * discount(r)
## Receive immediate profit if exercising
AOV[exercise] <- profit[t,exercise]

##Was the option exercised?
exercise_timing[exercise] <- t

###Re-iterate:
}
###End backwards induction

#Calculate project value: discount payoffs
option_values <- rep(0, G)

# American option value: discounting payoffs back to time zero, averaging over all paths.

exercised_paths <- !is.na(exercise_timing)
exercise_period <- exercise_timing[exercised_paths]
itm_paths <- which(exercised_paths)
option_values[exercised_paths] <- profit[nperiods * (itm_paths-1) + exercise_period] * discount(r, exercise_period)

# option value:
ROV <- mean(option_values)

if(!verbose)  return(ROV)

# Verbose Outputs:
exercise_time <- (exercise_period-1)*dt

return(list(
## Option value
`Value` = ROV,
## Option value standard error
`Standard Error` = sqrt(stats::var(option_values) / G),
## expected exercise time
`Expected Timing` = mean(exercise_time),
## exercise time standard error
`Expected Timing SE` = sqrt(stats::var(exercise_time) / G),
## exercise prob.
`Exercise Probability` = length(itm_paths) / G,
## cumulative exercise prob.
`Cumulative Exercise Probability` = cumsum(table(c(exercise_time, (0:(nperiods-1))*dt)) - 1) / G
))
}
```

## Try the LSMRealOptions package in your browser

Any scripts or data that you put into this service are public.

LSMRealOptions documentation built on June 26, 2021, 5:06 p.m.