Nothing
Option Pricing Functions to Accompany *Derivatives Markets*
In derivmkts: Functions and R Code to Accompany Derivatives Markets
rm(list=ls())
library(pander)
library(bookdown)
library(knitr)
library(ggplot2)
library(dplyr)
library(tidyr)
##homedir <- '/home/rmcd/tex/d67/Rtutorial/'
options(digits=4)
figsize <- 4.5
opts_chunk$set(size='footnotesize',
prompt=FALSE,
comment=NA
##,fig.align='center',
## fig.width = figsize,
## fig.height=figsize,
## out.width='3.75in'
)
opts_knit$set(#eval.after='fig.cap',
prompt=TRUE,
#renderer=renderer_latex(document=FALSE),
size='footnotesize')
curr <- function(amt) formatC(amt, format='f', digits=2)
library(derivmkts)
library(mnormt)
library(markdown)
opts_chunk$set(collapse=TRUE)
Introduction
This vignette is an overview to the functions in the derivmkts
package, which was conceived as a companion to my book
Derivatives Markets @mcdonald:derivs:13. The material
has an educational focus. There are other option pricing packages for
R, but this package has several distinguishing features:
- function names (mostly) correspond to those in \emph{Derivatives
Markets}.
- vectorized Greek calculations are convenient both for individual
options and for portfolios
- the
simprice function produces simulated asset prices
- the
quincunx function illustrates the workings of a
quincunx (Galton board).
- binomial functions include a plotting function that provides a
visual depiction of early exercise
Pricing functions and greeks
European Calls and Puts
Table \@ref(tab:bslist) lists the Black-Scholes
related functions in the package.^[See @black/scholes:73 and
@merton:73-bell.] The functions bscall, bsput, and bsopt
provide basic pricing of European calls and puts. There are also
options with binary payoffs: cash-or-nothing and asset-or-nothing
options. All of these functions are vectorized. The function bsopt
by default provides option greeks. Here are some examples:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0
bscall(s, k, v, r, tt, d)
bsput(s, c(95, 100, 105), v, r, tt, d)
bstbl <- data.frame(
Function = c('bscall', 'bsput', 'bsopt', 'assetcall', 'assetput',
'cashcall', 'cashput'),
Description = c('European call', 'European put', 'European call and put and associated Greeks: delta, gamma,
vega, theta, rho, psi, and elasticity', 'Asset-or-nothing call',
'Asset-or-nothing put', 'Cash-or-nothing call',
'Cash-or-nothing put'))
pander(bstbl,
caption = 'Black-Scholes related option pricing functions\\label{tab:bslist}',
split.cell = 60, justify = c('center', 'left'))
```{block2, include=FALSE}
\begin{table}[tp]
\centering
\begin{tabular}{cp{4in}}
Function& Description \ \hline
bscall & European call\
bsput & European put\
bsopt & European call and put and associated Greeks: delta, gamma,
vega, theta, rho, psi, and elasticity \
assetcall & Asset-or-nothing call\
assetput & Asset-or-nothing put\
cashcall & Cash-or-nothing call\
cashput & Cash-or-nothing put
\end{tabular}
\caption{Black-Scholes related option pricing functions}
\label{tab:bslist}
\end{table}
## Barrier Options
There are pricing functions for the following barrier
options:^[See @merton:73-bell, p. 175, for the
first derivation of a barrier option pricing formula and
@mcdonald:derivs:13, Chapter 14, for an overview.]
* down-and-in and down-and-out barrier binary options
* up-and-in and up-and-out barrier binary options
* more standard
down- and up- calls and puts, constructed using the barrier binary
options
Naming for the barrier options generally follows the convention
````
[u|d][i|o][call|put]
````
which means that the option is ``up'' or ``down'', ``in'' or ``out'', and a
call or put.^[This naming convention differs from that in
@mcdonald:derivs:13, in which names are `callupin`,
`callupout`, etc. Thus, I have made both names
available for these functions.] An up-and-in call, for example,
would be denoted by `uicall`. For binary options, we add the
underlying, which is either the asset or \$1: cash:
````
[asset|cash][u|d][i|o][call|put]
````
```r
H <- 115
bscall(s, c(80, 100, 120), v, r, tt, d)
## Up-and-in call
uicall(s, c(80, 100, 120), v, r, tt, d, H)
bsput(s, c(80, 100, 120), v, r, tt, d)
## Up-and-out put
uoput(s, c(80, 100, 120), v, r, tt, d, H)
Perpetual American Options
The functions callperpetual and putperetual
price infinitely-lived American
options.^[@merton:73-bell derived the price of a
perpetual American put.] The pricing formula assumes that all inputs
(risk-free rate, volatility, dividend yield) are fixed. This is of
course usual with the basic option pricing formulas, but it is more of
a conceptual stretch for an infinitely-lived option than for a 3-month
option.
In order for the option to have a determined value, the dividend yield
on the underlying asset must be positive if the option is a call. If
this is not true, the call is never exercised and the price is
undefined.^[A well-known result (see @merton:73-bell) is that
a standard American call is never exercised before expiration if the
dividend yield is zero and the interest rate is non-negative. A
perpetual call with $\delta=0$ and $r>0$ would thus never be
exercised. The limit of the option price as $\delta\to 0$ is $s$, so
in this case the function returns the stock price as the option
value.] Similarly, the risk-free rate must be positive if the option
is a put.
By default, the perpetual pricing formulas return the price. By
setting showbarrier=TRUE, the function returns both the option price
and the stock price at which the option is optimally exercised (the
``barrier''). Here are some examples:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.04
callperpetual(s, c(95, 100, 105), v, r, d)
callperpetual(s, c(95, 100, 105), v, r, d, showbarrier=TRUE)
Option Greeks
Options greeks are mathematical derivatives of the option price with
respect to inputs; see @mcdonald:derivs:13, Chapters 12 and 13, for a
discussion of the greeks for vanilla options. Greeks for vanilla and
barrier options can be computed using the greeks function, which is
a wrapper for any pricing function that returns the option price and
which uses the default naming of inputs.^[In this version of the
package, I have two alternative functions that return Greeks: a) The
bsopt function by default produces prices and Greeks for European
calls and puts, and b) The greeks2 function takes as arguments the
name of the pricing function and then inputs as a list. These may be
deprecated in the future. greeks2 is more cumbersome to use but may
be more robust. I welcome feedback on these functions and what you
find useful. ]
H <- 105
greeks(uicall(s, k, v, r, tt, d, H))
The value of this approach is that you can easily compute Greeks for
spreads and custom pricing functions. Here are two examples. First,
the value at time 0 of a prepaid contract that pays $S_{T}^{a}$ at
time $T$ is given by the powercontract() function:
powercontract <- function(s, v, r, tt, d, a) {
price <- exp(-r*tt)*s^a* exp((a*(r-d) + 1/2*a*(a-1)*v^2)*tt)
}
We can easily compute the Greeks for a power contract:
greeks(powercontract(s=40, v=.08, r=0.08, tt=0.25, d=0, a=2))
Second, consider a bull spread in which we buy a call with a strike of
$k_{1}$ and sell a call with a strike of $k_2$. We can create a
function that computes the
value of the spread, and then compute the greeks for the spread by
using this newly-created function together with greeks():
bullspread <- function(s, v, r, tt, d, k1, k2) {
bscall(s, k1, v, r, tt, d) - bscall(s, k2, v, r, tt, d)
}
greeks(bullspread(39:41, .3, .08, 1, 0, k1=40, k2=45))
\medskip
The Greeks function is vectorized, so you can create vectors of greek
values with a single call. This example plots, for a bull spread, the
gamma as a function of the stock price; see Figure
\@ref(fig:bullgamma).
\medskip
sseq <- seq(1, 100, by=0.5)
greeks(bullspread(sseq, .3, .08, 1, 0, k1=40, k2=45),
initcaps = TRUE, long = TRUE) %>%
filter(greek == 'Gamma' ) %>%
ggplot(aes(x = s, y = value)) + geom_line()
\medskip
This code produces the plots in Figure \@ref(fig:allgreeks):
\medskip
k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0
S <- seq(.5, 200, by=.5)
Call <- greeks(bscall(S, k, v, r, tt, d), long = TRUE)
Put <- greeks(bsput(S, k, v, r, tt, d), long = TRUE)
ggplot(rbind(Call, Put), aes(x = s, y = value, color = funcname )) +
geom_line() + facet_wrap(~ greek, scales = 'free_y') +
scale_color_discrete(name = 'Option', labels = c('Call','Put' )) +
scale_x_continuous('Stock', breaks =c(0, 100, 200) ) +
scale_y_continuous('Value')
Binomial Pricing of European and American Options
There are two functions related to binomial option
pricing:^[See @cox/ross/rubinstein:79,
@rendleman/bartter:79-jf, and @mcdonald:derivs:13, Chapter 11]
-
binomopt computes prices of American and European calls and
puts. The function has three optional parameters that control output:
-
returnparams=TRUE will return as a vector the option
pricing inputs, computed parameters, and risk-neutral probability.
-
returngreeks=TRUE will return as a vector the price,
delta, gamma, and theta at the initial node.
-
returntrees=TRUE will return as a list the price, greeks, the
full stock price tree, the exercise status (TRUE or FALSE) at
each node, and the replicating portfolio at each node.
-
binomplot displays the asset price
tree, the corresponding probability of being at each node, and
whether or not the option is exercised at each node. This
function is described in more detail in Section \@ref(sec:binomplot).
Here are examples of pricing, illustrating the default of just
returning the price, and the ability to return the price plus
parameters, as well as the price, the parameters, and various trees:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.03
binomopt(s, k, v, r, tt, d, nstep=3)
binomopt(s, k, v, r, tt, d, nstep=3, returnparams=TRUE)
binomopt(s, k, v, r, tt, d, nstep=3, putopt=TRUE)
binomopt(s, k, v, r, tt, d, nstep=3, returntrees=TRUE, putopt=TRUE)
Asian Options
There are analytical functions for valuing geometric Asian options and
Monte Carlo routines for valuing arithmetic Asian
options.^[See @kemna/vorst:90.] Be aware that the greeks()
function at this time will not work with the Monte Carlo valuation for
arithmetic Asian options. I plan to address this in a future
release.^[As the functions are currently written, each
invocation of the pricing function will start with a different random
number seed, resulting in price variation that is due solely to random
variation. Moreover, random number generation changes the random
number seed globally. In a future release I hope to address this by
saving and restoring the seed within the greeks function. For the
curious, a Stackoverflow
post
discusses this issue.]
Geometric Asian Options
Geometric Asian options can be valued using the Black-Scholes formulas
for vanilla calls and puts, with modified inputs. The functions return
both call and put prices with a named vector:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.03; m <- 3
geomavgpricecall(s, 98:102, v, r, tt, d, m)
geomavgpricecall(s, 98:102, v, r, tt, d, m, cont=TRUE)
geomavgstrikecall(s, k, v, r, tt, d, m)
Arithmetic Asian Options
Monte Carlo valuation is used to price arithmetic Asian options. For
efficiency, the function arithasianmc returns call and put
prices for average price and average strike options. By default the
number of simulations is 1000. Optionally the function returns the
standard deviation of each estimate
arithasianmc(s, k, v, r, tt, d, 3, numsim=5000, printsds=TRUE)
The function arithavgpricecv uses the control variate
method to reduce the variance in the simulation. At the moment this
function prices only calls, and returns both the price and the
regression coefficient used in the control variate correction:
arithavgpricecv(s, k, v, r, tt, d, 3, numsim=5000)
Compound Options
A compound option is an option where the underlying asset is an
option.^[See @geske:79 and @mcdonald:derivs:13, Chapter 14.]
The terminology associated with compound options can be confusing, so
it may be easiest to start with an example.
Figure \@ref(fig:compoundopt) is a timeline for a compound
option that is an option to buy a put. The compound option expires at
$t_{1}$ and the put expires at $t_{2}$. The owner of the compound
option only acquires the put if at time $t_{1}$ it is worth at least
$k_{co}$, and only exercises the put if at time $t_{2}$ the stock
price is less than $k_{uo}$.
compound.caption <- 'The timeline for a compound option: a call to buy a put. The compound option expires at time $t_{1}$ and the underlying asset is a put option that expires at time $t_{2}$. At time $t_{1}$, the owner decides whether to pay $k_{co}$ to buy a put option which has time to expiration $t_{2} - t_{1}$. At time $t_{2}$ the owner decides whether to exercise the put, selling the stock for the strike price of $k_{uo}$.'
# \label{fig:compoundopt}
#\end{figure}
```r$ and the underlying asset is a put option that expires at time $t_{2}$. At time $t_{1}$, the owner decides whether to pay $k_{co}$ to buy a put option which has time to expiration $t_{2} - t_{1}$. At time $t_{2}$ the owner decides whether to exercise the put, selling the stock for the strike price of $k_{uo}$.', echo=FALSE}
\begin{tikzpicture}
\newcommand{\start}{0}
\newcommand{\finish}{8}
\newcommand{\midpt}{(\finish/2-\start/2}
\newcommand{\tickheight}{0.3}
\draw (\start,0) -- (\finish,0);
\draw (0,\tickheight) -- (0, -\tickheight)
node[below] {\begin{tabular}{c} Time 0 \ Buy compound \option
\end{tabular}};
\draw (\midpt,\tickheight) -- (\midpt, -\tickheight)
node[below]
{\begin{tabular}{c}
Time $t_{1}$ \
Compound exercise \
decision: Pay $k_{co}$ \
to buy put?
\end{tabular}
};
\draw (\finish, \tickheight) -- (\finish, -\tickheight)
node[below] {
\begin{tabular}{c}
Time $t_{2}$ \
Put exercise decision\
Sell stock for $k_{uo}$?
\end{tabular}
};
\end{tikzpicture}
### Definition of a Compound Option
Based on the example, you can see that there are three prices
associated with a compound option:
* The price of an underlying asset.
* The price of the underlying option, which is an option to buy or
sell the underlying asset (we will refer to this as the price of
the underlying option)
* The price of the compound option, which gives us the right to
buy or sell the underlying option
The definition of a compound option therefore requires that we specify
* whether the underlying option is a put or a call
* whether the compound option is a put or a call
* the strike price at which you can exercise the underlying option
($k_{uo}$)
* the strike price at which you can exercise the compound option ($k_{co}$)
* the date at which you can exercise the compound option (first
exercise date, $t_{1}$)
* the date at which you can exercise the underlying option expires, $t_{2}>t_{1}$.
Given these possibilities, you can have a call on a call, a put on a
call, a call on a put, and a put on a put. The valuation procedure
require calculating the underlying asset price at which you are
indifferent about acquiring the underlying option.
The price calculation requires computing the stock price above or
below which you would optimally exercise the option at time
$t_{1}$.
### Examples
As an example, consider the following inputs for a call option to buy
a call option:
```r
s <- 100; kuo <- 95; v <- 0.30; r <- 0.08; t1 <- 0.50; t2 <- 0.75; d <- 0
kco <- 3.50
calloncall(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
With these parameters, after 6 months ($t_{1}=0.5$), the compound
option buyer decides whether to pay \$r formatC(kco, format='f',
digits=2) to acquire a 3-month call on the underlying asset. (The
volatility of the underlying asset is r v.) It will be
worthwhile to pay the compound strike, \$r curr(kco), as long as
the underlying asset price exceeds r calloncall(s, kuo, kco, v,
r, t1, t2, d, returnscritical=TRUE)['scritical'].
Similarly, there is a put on the call, and a call and put on the
corresponding put:
putoncall(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
callonput(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
putonput(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
Jumps and Stochastic Volatility
\label{sec:jumps}
The mertonjump function returns call and put prices for a
stock that can jump discretely.^[See @merton/jfe/1976.]
A poisson process controls the
occurrence of a jump and the size of the jump is lognormally
distributed. The parameter lambda is the mean number of
jumps per year, the parameter alphaj is the log of the
expected jump, and sigmaj is the standard deviation of the
log of the jump. The jump amount is thus drawn from the distribution
\begin{equation}
Y \sim \mathcal{N}(\alpha_{J} - 0.5\sigma^{2}{J}, \sigma{J}^{2} )
\end{equation}
mertonjump(s, k, v, r, tt, d, lambda=0.5, alphaj=-0.2, vj=0.3)
c(bscall(s, k, v, r, tt, d), bsput(s, k, v, r, tt, d))
Bonds
The simple bond functions provided in this version compute the present
value of cash flows (bondpv), the IRR of the bond (bondyield),
Macaulay duration (duration), and convexity (convexity).
coupon <- 8; mat <- 20; yield <- 0.06; principal <- 100;
modified <- FALSE; freq <- 2
price <- bondpv(coupon, mat, yield, principal, freq)
price
bondyield(price, coupon, mat, principal, freq)
duration(price, coupon, mat, principal, freq, modified)
convexity(price, coupon, mat, principal, freq)
Monte Carlo simulation of prices
The function simprice provides a flexible way to generate prices
that can be used for Monte Carlo valuation or just to illustrate
sample paths.^[For more information on Monte Carlo see
@mcdonald:derivs:13, Chapter 19, and @glasserman:mc:04, especially
Chapter 4.] This function implements "naive" Monte Carlo, not using
any variance reduction techniques. When supplied with a covariance
matrix, simprice produces simulations of multiple assets with
correlated returns. There are default values for all inputs.
Other option pricing functions in this package assume that you specify
volatility as the annualized return standard deviation. With a scalar
volatility input, the simprice function also interprets the value as
a standard deviation. A matrix, however, is interpreted as a
covariance matrix, which means that the individual stock volatilities
are square root of the diagonal elements. Because scalar and matrix
inputs are treated differently, there is an option To change the
behavior for a scalar input. To interpret the value as a variance,
specify scalar_v_is_stddev = FALSE.
Long vs wide output
The function simprice can produce either long or wide output. Here
are examples of each, using the default parameters.^[The simprice
function save and restores the random number seed, so repeated
invocations without intermediate seed-changing actions will produce
the same result.]
args(simprice)
simprice(long = TRUE)
simprice(long = FALSE)
Simulated price paths
A simple example is to generate 5 random sample paths of daily stock
prices for a year; see Figure \@ref(fig:fivepaths):
s0 <- 100; v <- 0.3; r <- 0.10; d <- 0; tt <- 1
trials <- 5; periods <- 365; set.seed(1)
s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials,
periods = periods, jump = FALSE, long = TRUE)
ggplot(s, aes(x = period, y = price, color = trial, group = trial)) +
geom_line()
We can do the same exercise adding jumps; see Figure \@ref(fig:fivejumpers).
s0 <- 100; v <- 0.3; r <- 0.10; d <- 0; tt <- 1
trials <- 5; periods <- 365; jump <- TRUE; lambda <- 2;
alphaj <- -0.1; vj <- 0.2; set.seed(1)
s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials,
periods = periods, jump = jump, alphaj = alphaj,
lambda = lambda, vj = vj, long = TRUE)
ggplot(s, aes(x = period, y = price, color = trial, group = trial)) +
geom_line()
\newpage
Multiple correlated stocks
When you supply an $n\times n$ covariance matrix as the volatility input,
simprice produces simulations for $n$ stocks using the specified
covariances. This makes it convenient to use an esimated covariance
matrix to produce a simulation.
Negatively correlated assets
To illustrate the use of a covariance matrix, Figure
\@ref(fig:negcoorsim) plots two stocks which are highly negatively
correlated.
set.seed(1)
s0 <- 100; r <- 0.08; tt <- 1; d <- 0; jump <- FALSE
trials <- 1; periods <- 52;
v <- .3^2*matrix(c(1, -.99, -.99, 1), nrow = 2)
s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials,
periods = periods, jump = jump, long = TRUE)
ggplot(s, aes(x = period, y = price, group = asset, color = asset)) +
geom_line()
Three correlated assets
Here is an example using a three-asset covariance matrix, showing that
the simulated series have the correct variance structure. We have
three correlated stocks with standard deviations of 60%, 20%, and
45%., with the covariance matrix given by v. We can compare the
estimated covariance matrix with the input,v:
set.seed(1)
tt <- 2; periods <- tt*365
vc <- vols <- diag(3)
diag(vols) <- c(.6, .2, .45) ## volatilities
corrs <- c(.4, -.3, .25)
vc[lower.tri(vc)] <- corrs ## correlations
vc <- t(vc) ## lower triangular becomes upper triangular
vc[lower.tri(vc)] <- corrs
v <- vols %*% vc %*% vols
v
s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = 1,
periods = periods, jump = FALSE, long = TRUE)
threestocks <- s %>%
filter(trial == 1) %>%
group_by(asset) %>%
mutate(ret = log(price/lag(price)),
row = row_number()) %>%
select(asset, period, ret) %>%
spread(key = asset, value = ret )
var(threestocks[2:4], na.rm = TRUE)*365
Functions with Graphical Output
Several functions provide visual illustrations of some aspects of the
material.
Quincunx or Galton Board
The quincunx is a physical device the illustrates the central limit
theorem. A ball rolls down a pegboard and strikes a peg, falling
randomly either to the left or right. As it continues down the board
it continues to strike a series of pegs, randomly falling left or
right at each. The balls
collect in bins and create an approximate normal distribution.
The quincunx function allows the user to simulate a quincunx,
observing the path of each ball and watching the height of each bin as
the balls accumulate. More interestingly, the quincunx function
permits altering the probability that the ball will fall to the
right.
Figure \@ref(fig:quincunx) illustrates the function after dropping 200
balls down 20 levels of pegs with a 70\% probability that each ball
will fall right:
par(mar=c(2,2,2,2))
quincunx(n=20, numballs=200, delay=0, probright=0.7)
Plotting the Solution to the Binomial Pricing Model {#sec:binomplot}
\label{sec:binomplot}
The binomplot function calls binomopt to compute the option price
and the various trees, which it then uses in plotting:
The first plot, figure \@ref(fig:binomplot1), is basic:
s0 <- 100; k <- 100; v <- 0.3; r <- 0.08; tt <- 2; d <- 0
binomplot(s0, k, v, r, tt, d, nstep=6, american=TRUE, putopt=TRUE)
The second plot, figure \@ref(fig:binomplot2), adds a display of stock
prices and arrows connecting the nodes.
binomplot(s0, k, v, r, tt, d, nstep=6, american=TRUE, putopt=TRUE,
plotvalues=TRUE, plotarrows=TRUE)
As a final example, consider an American call when the dividend yield
is positive and nstep has a larger value. Figure
\@ref(fig:binomplot3) shows the plot, with early exercise evident.
d <- 0.06
binomplot(s0, k, v, r, tt, d, nstep=40, american=TRUE)
The large value of nstep creates a high maximum terminal stock
price, which makes details hard to discern in the boundary region
where exercise first occurrs. We can zoom in on that region by
selecting values for ylimval; the result is in Figure \@ref(fig:binomplot4).
d <- 0.06
binomplot(s0, k, v, r, tt, d, nstep=40, american=TRUE, ylimval=c(75, 225))
References
Try the derivmkts package in your browser
Any scripts or data that you put into this service are public.
derivmkts documentation built on April 11, 2022, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
rm(list=ls()) library(pander) library(bookdown) library(knitr) library(ggplot2) library(dplyr) library(tidyr) ##homedir <- '/home/rmcd/tex/d67/Rtutorial/' options(digits=4) figsize <- 4.5 opts_chunk$set(size='footnotesize', prompt=FALSE, comment=NA ##,fig.align='center', ## fig.width = figsize, ## fig.height=figsize, ## out.width='3.75in' ) opts_knit$set(#eval.after='fig.cap', prompt=TRUE, #renderer=renderer_latex(document=FALSE), size='footnotesize') curr <- function(amt) formatC(amt, format='f', digits=2)
library(derivmkts) library(mnormt) library(markdown) opts_chunk$set(collapse=TRUE)
Introduction
This vignette is an overview to the functions in the derivmkts package, which was conceived as a companion to my book Derivatives Markets @mcdonald:derivs:13. The material has an educational focus. There are other option pricing packages for R, but this package has several distinguishing features:
- function names (mostly) correspond to those in \emph{Derivatives Markets}.
- vectorized Greek calculations are convenient both for individual options and for portfolios
- the
simpricefunction produces simulated asset prices - the
quincunxfunction illustrates the workings of a quincunx (Galton board). - binomial functions include a plotting function that provides a visual depiction of early exercise
Pricing functions and greeks
European Calls and Puts
Table \@ref(tab:bslist) lists the Black-Scholes
related functions in the package.^[See @black/scholes:73 and
@merton:73-bell.] The functions bscall, bsput, and bsopt
provide basic pricing of European calls and puts. There are also
options with binary payoffs: cash-or-nothing and asset-or-nothing
options. All of these functions are vectorized. The function bsopt
by default provides option greeks. Here are some examples:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0 bscall(s, k, v, r, tt, d) bsput(s, c(95, 100, 105), v, r, tt, d)
bstbl <- data.frame( Function = c('bscall', 'bsput', 'bsopt', 'assetcall', 'assetput', 'cashcall', 'cashput'), Description = c('European call', 'European put', 'European call and put and associated Greeks: delta, gamma, vega, theta, rho, psi, and elasticity', 'Asset-or-nothing call', 'Asset-or-nothing put', 'Cash-or-nothing call', 'Cash-or-nothing put')) pander(bstbl, caption = 'Black-Scholes related option pricing functions\\label{tab:bslist}', split.cell = 60, justify = c('center', 'left'))
```{block2, include=FALSE} \begin{table}[tp] \centering \begin{tabular}{cp{4in}} Function& Description \ \hline bscall & European call\ bsput & European put\ bsopt & European call and put and associated Greeks: delta, gamma, vega, theta, rho, psi, and elasticity \ assetcall & Asset-or-nothing call\ assetput & Asset-or-nothing put\ cashcall & Cash-or-nothing call\ cashput & Cash-or-nothing put \end{tabular} \caption{Black-Scholes related option pricing functions} \label{tab:bslist} \end{table}
## Barrier Options There are pricing functions for the following barrier options:^[See @merton:73-bell, p. 175, for the first derivation of a barrier option pricing formula and @mcdonald:derivs:13, Chapter 14, for an overview.] * down-and-in and down-and-out barrier binary options * up-and-in and up-and-out barrier binary options * more standard down- and up- calls and puts, constructed using the barrier binary options Naming for the barrier options generally follows the convention ```` [u|d][i|o][call|put] ```` which means that the option is ``up'' or ``down'', ``in'' or ``out'', and a call or put.^[This naming convention differs from that in @mcdonald:derivs:13, in which names are `callupin`, `callupout`, etc. Thus, I have made both names available for these functions.] An up-and-in call, for example, would be denoted by `uicall`. For binary options, we add the underlying, which is either the asset or \$1: cash: ```` [asset|cash][u|d][i|o][call|put] ```` ```r H <- 115 bscall(s, c(80, 100, 120), v, r, tt, d) ## Up-and-in call uicall(s, c(80, 100, 120), v, r, tt, d, H) bsput(s, c(80, 100, 120), v, r, tt, d) ## Up-and-out put uoput(s, c(80, 100, 120), v, r, tt, d, H)
Perpetual American Options
The functions callperpetual and putperetual
price infinitely-lived American
options.^[@merton:73-bell derived the price of a
perpetual American put.] The pricing formula assumes that all inputs
(risk-free rate, volatility, dividend yield) are fixed. This is of
course usual with the basic option pricing formulas, but it is more of
a conceptual stretch for an infinitely-lived option than for a 3-month
option.
In order for the option to have a determined value, the dividend yield on the underlying asset must be positive if the option is a call. If this is not true, the call is never exercised and the price is undefined.^[A well-known result (see @merton:73-bell) is that a standard American call is never exercised before expiration if the dividend yield is zero and the interest rate is non-negative. A perpetual call with $\delta=0$ and $r>0$ would thus never be exercised. The limit of the option price as $\delta\to 0$ is $s$, so in this case the function returns the stock price as the option value.] Similarly, the risk-free rate must be positive if the option is a put.
By default, the perpetual pricing formulas return the price. By
setting showbarrier=TRUE, the function returns both the option price
and the stock price at which the option is optimally exercised (the
``barrier''). Here are some examples:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.04 callperpetual(s, c(95, 100, 105), v, r, d) callperpetual(s, c(95, 100, 105), v, r, d, showbarrier=TRUE)
Option Greeks
Options greeks are mathematical derivatives of the option price with
respect to inputs; see @mcdonald:derivs:13, Chapters 12 and 13, for a
discussion of the greeks for vanilla options. Greeks for vanilla and
barrier options can be computed using the greeks function, which is
a wrapper for any pricing function that returns the option price and
which uses the default naming of inputs.^[In this version of the
package, I have two alternative functions that return Greeks: a) The
bsopt function by default produces prices and Greeks for European
calls and puts, and b) The greeks2 function takes as arguments the
name of the pricing function and then inputs as a list. These may be
deprecated in the future. greeks2 is more cumbersome to use but may
be more robust. I welcome feedback on these functions and what you
find useful. ]
H <- 105 greeks(uicall(s, k, v, r, tt, d, H))
The value of this approach is that you can easily compute Greeks for
spreads and custom pricing functions. Here are two examples. First,
the value at time 0 of a prepaid contract that pays $S_{T}^{a}$ at
time $T$ is given by the powercontract() function:
powercontract <- function(s, v, r, tt, d, a) { price <- exp(-r*tt)*s^a* exp((a*(r-d) + 1/2*a*(a-1)*v^2)*tt) }
We can easily compute the Greeks for a power contract:
greeks(powercontract(s=40, v=.08, r=0.08, tt=0.25, d=0, a=2))
Second, consider a bull spread in which we buy a call with a strike of
$k_{1}$ and sell a call with a strike of $k_2$. We can create a
function that computes the
value of the spread, and then compute the greeks for the spread by
using this newly-created function together with greeks():
bullspread <- function(s, v, r, tt, d, k1, k2) { bscall(s, k1, v, r, tt, d) - bscall(s, k2, v, r, tt, d) } greeks(bullspread(39:41, .3, .08, 1, 0, k1=40, k2=45))
\medskip
The Greeks function is vectorized, so you can create vectors of greek values with a single call. This example plots, for a bull spread, the gamma as a function of the stock price; see Figure \@ref(fig:bullgamma).
\medskip
sseq <- seq(1, 100, by=0.5) greeks(bullspread(sseq, .3, .08, 1, 0, k1=40, k2=45), initcaps = TRUE, long = TRUE) %>% filter(greek == 'Gamma' ) %>% ggplot(aes(x = s, y = value)) + geom_line()
\medskip
This code produces the plots in Figure \@ref(fig:allgreeks):
\medskip
k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0 S <- seq(.5, 200, by=.5) Call <- greeks(bscall(S, k, v, r, tt, d), long = TRUE) Put <- greeks(bsput(S, k, v, r, tt, d), long = TRUE) ggplot(rbind(Call, Put), aes(x = s, y = value, color = funcname )) + geom_line() + facet_wrap(~ greek, scales = 'free_y') + scale_color_discrete(name = 'Option', labels = c('Call','Put' )) + scale_x_continuous('Stock', breaks =c(0, 100, 200) ) + scale_y_continuous('Value')
Binomial Pricing of European and American Options
There are two functions related to binomial option pricing:^[See @cox/ross/rubinstein:79, @rendleman/bartter:79-jf, and @mcdonald:derivs:13, Chapter 11]
-
binomopt computes prices of American and European calls and puts. The function has three optional parameters that control output:
-
returnparams=TRUEwill return as a vector the option pricing inputs, computed parameters, and risk-neutral probability. -
returngreeks=TRUEwill return as a vector the price, delta, gamma, and theta at the initial node. -
returntrees=TRUEwill return as a list the price, greeks, the full stock price tree, the exercise status (TRUEorFALSE) at each node, and the replicating portfolio at each node. -
binomplot displays the asset price tree, the corresponding probability of being at each node, and whether or not the option is exercised at each node. This function is described in more detail in Section \@ref(sec:binomplot).
Here are examples of pricing, illustrating the default of just returning the price, and the ability to return the price plus parameters, as well as the price, the parameters, and various trees:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.03 binomopt(s, k, v, r, tt, d, nstep=3) binomopt(s, k, v, r, tt, d, nstep=3, returnparams=TRUE) binomopt(s, k, v, r, tt, d, nstep=3, putopt=TRUE) binomopt(s, k, v, r, tt, d, nstep=3, returntrees=TRUE, putopt=TRUE)
Asian Options
There are analytical functions for valuing geometric Asian options and
Monte Carlo routines for valuing arithmetic Asian
options.^[See @kemna/vorst:90.] Be aware that the greeks()
function at this time will not work with the Monte Carlo valuation for
arithmetic Asian options. I plan to address this in a future
release.^[As the functions are currently written, each
invocation of the pricing function will start with a different random
number seed, resulting in price variation that is due solely to random
variation. Moreover, random number generation changes the random
number seed globally. In a future release I hope to address this by
saving and restoring the seed within the greeks function. For the
curious, a Stackoverflow
post
discusses this issue.]
Geometric Asian Options
Geometric Asian options can be valued using the Black-Scholes formulas for vanilla calls and puts, with modified inputs. The functions return both call and put prices with a named vector:
s <- 100; k <- 100; r <- 0.08; v <- 0.30; tt <- 2; d <- 0.03; m <- 3 geomavgpricecall(s, 98:102, v, r, tt, d, m) geomavgpricecall(s, 98:102, v, r, tt, d, m, cont=TRUE) geomavgstrikecall(s, k, v, r, tt, d, m)
Arithmetic Asian Options
Monte Carlo valuation is used to price arithmetic Asian options. For
efficiency, the function arithasianmc returns call and put
prices for average price and average strike options. By default the
number of simulations is 1000. Optionally the function returns the
standard deviation of each estimate
arithasianmc(s, k, v, r, tt, d, 3, numsim=5000, printsds=TRUE)
The function arithavgpricecv uses the control variate
method to reduce the variance in the simulation. At the moment this
function prices only calls, and returns both the price and the
regression coefficient used in the control variate correction:
arithavgpricecv(s, k, v, r, tt, d, 3, numsim=5000)
Compound Options
A compound option is an option where the underlying asset is an option.^[See @geske:79 and @mcdonald:derivs:13, Chapter 14.] The terminology associated with compound options can be confusing, so it may be easiest to start with an example.
Figure \@ref(fig:compoundopt) is a timeline for a compound option that is an option to buy a put. The compound option expires at $t_{1}$ and the put expires at $t_{2}$. The owner of the compound option only acquires the put if at time $t_{1}$ it is worth at least $k_{co}$, and only exercises the put if at time $t_{2}$ the stock price is less than $k_{uo}$.
compound.caption <- 'The timeline for a compound option: a call to buy a put. The compound option expires at time $t_{1}$ and the underlying asset is a put option that expires at time $t_{2}$. At time $t_{1}$, the owner decides whether to pay $k_{co}$ to buy a put option which has time to expiration $t_{2} - t_{1}$. At time $t_{2}$ the owner decides whether to exercise the put, selling the stock for the strike price of $k_{uo}$.' # \label{fig:compoundopt} #\end{figure}
```r$ and the underlying asset is a put option that expires at time $t_{2}$. At time $t_{1}$, the owner decides whether to pay $k_{co}$ to buy a put option which has time to expiration $t_{2} - t_{1}$. At time $t_{2}$ the owner decides whether to exercise the put, selling the stock for the strike price of $k_{uo}$.', echo=FALSE} \begin{tikzpicture} \newcommand{\start}{0} \newcommand{\finish}{8} \newcommand{\midpt}{(\finish/2-\start/2} \newcommand{\tickheight}{0.3} \draw (\start,0) -- (\finish,0); \draw (0,\tickheight) -- (0, -\tickheight) node[below] {\begin{tabular}{c} Time 0 \ Buy compound \option \end{tabular}}; \draw (\midpt,\tickheight) -- (\midpt, -\tickheight) node[below] {\begin{tabular}{c} Time $t_{1}$ \ Compound exercise \ decision: Pay $k_{co}$ \ to buy put? \end{tabular} }; \draw (\finish, \tickheight) -- (\finish, -\tickheight) node[below] { \begin{tabular}{c} Time $t_{2}$ \ Put exercise decision\ Sell stock for $k_{uo}$? \end{tabular} }; \end{tikzpicture}
### Definition of a Compound Option Based on the example, you can see that there are three prices associated with a compound option: * The price of an underlying asset. * The price of the underlying option, which is an option to buy or sell the underlying asset (we will refer to this as the price of the underlying option) * The price of the compound option, which gives us the right to buy or sell the underlying option The definition of a compound option therefore requires that we specify * whether the underlying option is a put or a call * whether the compound option is a put or a call * the strike price at which you can exercise the underlying option ($k_{uo}$) * the strike price at which you can exercise the compound option ($k_{co}$) * the date at which you can exercise the compound option (first exercise date, $t_{1}$) * the date at which you can exercise the underlying option expires, $t_{2}>t_{1}$. Given these possibilities, you can have a call on a call, a put on a call, a call on a put, and a put on a put. The valuation procedure require calculating the underlying asset price at which you are indifferent about acquiring the underlying option. The price calculation requires computing the stock price above or below which you would optimally exercise the option at time $t_{1}$. ### Examples As an example, consider the following inputs for a call option to buy a call option: ```r s <- 100; kuo <- 95; v <- 0.30; r <- 0.08; t1 <- 0.50; t2 <- 0.75; d <- 0 kco <- 3.50 calloncall(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
With these parameters, after 6 months ($t_{1}=0.5$), the compound
option buyer decides whether to pay \$r formatC(kco, format='f',
digits=2) to acquire a 3-month call on the underlying asset. (The
volatility of the underlying asset is r v.) It will be
worthwhile to pay the compound strike, \$r curr(kco), as long as
the underlying asset price exceeds r calloncall(s, kuo, kco, v,
r, t1, t2, d, returnscritical=TRUE)['scritical'].
Similarly, there is a put on the call, and a call and put on the corresponding put:
putoncall(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE) callonput(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE) putonput(s, kuo, kco, v, r, t1, t2, d, returnscritical=TRUE)
Jumps and Stochastic Volatility
\label{sec:jumps}
The mertonjump function returns call and put prices for a
stock that can jump discretely.^[See @merton/jfe/1976.]
A poisson process controls the
occurrence of a jump and the size of the jump is lognormally
distributed. The parameter lambda is the mean number of
jumps per year, the parameter alphaj is the log of the
expected jump, and sigmaj is the standard deviation of the
log of the jump. The jump amount is thus drawn from the distribution
\begin{equation}
Y \sim \mathcal{N}(\alpha_{J} - 0.5\sigma^{2}{J}, \sigma{J}^{2} )
\end{equation}
mertonjump(s, k, v, r, tt, d, lambda=0.5, alphaj=-0.2, vj=0.3) c(bscall(s, k, v, r, tt, d), bsput(s, k, v, r, tt, d))
Bonds
The simple bond functions provided in this version compute the present
value of cash flows (bondpv), the IRR of the bond (bondyield),
Macaulay duration (duration), and convexity (convexity).
coupon <- 8; mat <- 20; yield <- 0.06; principal <- 100; modified <- FALSE; freq <- 2 price <- bondpv(coupon, mat, yield, principal, freq) price bondyield(price, coupon, mat, principal, freq) duration(price, coupon, mat, principal, freq, modified) convexity(price, coupon, mat, principal, freq)
Monte Carlo simulation of prices
The function simprice provides a flexible way to generate prices
that can be used for Monte Carlo valuation or just to illustrate
sample paths.^[For more information on Monte Carlo see
@mcdonald:derivs:13, Chapter 19, and @glasserman:mc:04, especially
Chapter 4.] This function implements "naive" Monte Carlo, not using
any variance reduction techniques. When supplied with a covariance
matrix, simprice produces simulations of multiple assets with
correlated returns. There are default values for all inputs.
Other option pricing functions in this package assume that you specify
volatility as the annualized return standard deviation. With a scalar
volatility input, the simprice function also interprets the value as
a standard deviation. A matrix, however, is interpreted as a
covariance matrix, which means that the individual stock volatilities
are square root of the diagonal elements. Because scalar and matrix
inputs are treated differently, there is an option To change the
behavior for a scalar input. To interpret the value as a variance,
specify scalar_v_is_stddev = FALSE.
Long vs wide output
The function simprice can produce either long or wide output. Here
are examples of each, using the default parameters.^[The simprice
function save and restores the random number seed, so repeated
invocations without intermediate seed-changing actions will produce
the same result.]
args(simprice) simprice(long = TRUE) simprice(long = FALSE)
Simulated price paths
A simple example is to generate 5 random sample paths of daily stock prices for a year; see Figure \@ref(fig:fivepaths):
s0 <- 100; v <- 0.3; r <- 0.10; d <- 0; tt <- 1 trials <- 5; periods <- 365; set.seed(1) s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials, periods = periods, jump = FALSE, long = TRUE) ggplot(s, aes(x = period, y = price, color = trial, group = trial)) + geom_line()
We can do the same exercise adding jumps; see Figure \@ref(fig:fivejumpers).
s0 <- 100; v <- 0.3; r <- 0.10; d <- 0; tt <- 1 trials <- 5; periods <- 365; jump <- TRUE; lambda <- 2; alphaj <- -0.1; vj <- 0.2; set.seed(1) s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials, periods = periods, jump = jump, alphaj = alphaj, lambda = lambda, vj = vj, long = TRUE) ggplot(s, aes(x = period, y = price, color = trial, group = trial)) + geom_line()
\newpage
Multiple correlated stocks
When you supply an $n\times n$ covariance matrix as the volatility input,
simprice produces simulations for $n$ stocks using the specified
covariances. This makes it convenient to use an esimated covariance
matrix to produce a simulation.
Negatively correlated assets
To illustrate the use of a covariance matrix, Figure \@ref(fig:negcoorsim) plots two stocks which are highly negatively correlated.
set.seed(1) s0 <- 100; r <- 0.08; tt <- 1; d <- 0; jump <- FALSE trials <- 1; periods <- 52; v <- .3^2*matrix(c(1, -.99, -.99, 1), nrow = 2) s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = trials, periods = periods, jump = jump, long = TRUE) ggplot(s, aes(x = period, y = price, group = asset, color = asset)) + geom_line()
Three correlated assets
Here is an example using a three-asset covariance matrix, showing that
the simulated series have the correct variance structure. We have
three correlated stocks with standard deviations of 60%, 20%, and
45%., with the covariance matrix given by v. We can compare the
estimated covariance matrix with the input,v:
set.seed(1) tt <- 2; periods <- tt*365 vc <- vols <- diag(3) diag(vols) <- c(.6, .2, .45) ## volatilities corrs <- c(.4, -.3, .25) vc[lower.tri(vc)] <- corrs ## correlations vc <- t(vc) ## lower triangular becomes upper triangular vc[lower.tri(vc)] <- corrs v <- vols %*% vc %*% vols v s <- simprice(s0 = s0, v = v, r = r, tt = tt, d = d, trials = 1, periods = periods, jump = FALSE, long = TRUE) threestocks <- s %>% filter(trial == 1) %>% group_by(asset) %>% mutate(ret = log(price/lag(price)), row = row_number()) %>% select(asset, period, ret) %>% spread(key = asset, value = ret ) var(threestocks[2:4], na.rm = TRUE)*365
Functions with Graphical Output
Several functions provide visual illustrations of some aspects of the material.
Quincunx or Galton Board
The quincunx is a physical device the illustrates the central limit theorem. A ball rolls down a pegboard and strikes a peg, falling randomly either to the left or right. As it continues down the board it continues to strike a series of pegs, randomly falling left or right at each. The balls collect in bins and create an approximate normal distribution.
The quincunx function allows the user to simulate a quincunx, observing the path of each ball and watching the height of each bin as the balls accumulate. More interestingly, the quincunx function permits altering the probability that the ball will fall to the right.
Figure \@ref(fig:quincunx) illustrates the function after dropping 200 balls down 20 levels of pegs with a 70\% probability that each ball will fall right:
par(mar=c(2,2,2,2)) quincunx(n=20, numballs=200, delay=0, probright=0.7)
Plotting the Solution to the Binomial Pricing Model {#sec:binomplot}
\label{sec:binomplot}
The binomplot function calls binomopt to compute the option price
and the various trees, which it then uses in plotting:
The first plot, figure \@ref(fig:binomplot1), is basic:
s0 <- 100; k <- 100; v <- 0.3; r <- 0.08; tt <- 2; d <- 0 binomplot(s0, k, v, r, tt, d, nstep=6, american=TRUE, putopt=TRUE)
The second plot, figure \@ref(fig:binomplot2), adds a display of stock prices and arrows connecting the nodes.
binomplot(s0, k, v, r, tt, d, nstep=6, american=TRUE, putopt=TRUE, plotvalues=TRUE, plotarrows=TRUE)
As a final example, consider an American call when the dividend yield
is positive and nstep has a larger value. Figure
\@ref(fig:binomplot3) shows the plot, with early exercise evident.
d <- 0.06 binomplot(s0, k, v, r, tt, d, nstep=40, american=TRUE)
The large value of nstep creates a high maximum terminal stock
price, which makes details hard to discern in the boundary region
where exercise first occurrs. We can zoom in on that region by
selecting values for ylimval; the result is in Figure \@ref(fig:binomplot4).
d <- 0.06 binomplot(s0, k, v, r, tt, d, nstep=40, american=TRUE, ylimval=c(75, 225))
References
Try the derivmkts package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.