Type:  Package 
Title:  Numerical Methods and Optimization in Finance 
Version:  2.80 
Date:  20231020 
Maintainer:  Enrico Schumann <es@enricoschumann.net> 
Depends:  R (≥ 3.5) 
Imports:  grDevices, graphics, parallel, stats, utils 
Suggests:  MASS, PMwR, RUnit, Rglpk, datetimeutils, openxlsx, quadprog, readxl, tinytest 
Description:  Functions, examples and data from the first and the second edition of "Numerical Methods and Optimization in Finance" by M. Gilli, D. Maringer and E. Schumann (2019, ISBN:9780128150658). The package provides implementations of optimisation heuristics (Differential Evolution, Genetic Algorithms, Particle Swarm Optimisation, Simulated Annealing and Threshold Accepting), and other optimisation tools, such as grid search and greedy search. There are also functions for the valuation of financial instruments such as bonds and options, for portfolio selection and functions that help with stochastic simulations. 
License:  GPL3 
URL:  http://enricoschumann.net/NMOF.htm , https://gitlab.com/NMOF , https://git.sr.ht/~enricoschumann/NMOF , https://github.com/enricoschumann/NMOF 
LazyLoad:  yes 
LazyData:  yes 
ByteCompile:  yes 
Classification/JEL:  C61, C63 
NeedsCompilation:  no 
Packaged:  20231020 07:49:45 UTC; es19 
Author:  Enrico Schumann [aut, cre] 
Repository:  CRAN 
Date/Publication:  20231020 08:20:02 UTC 
Built:  R 4.5.0; ; 20240430 19:21:07 UTC; unix 
Numerical Methods and Optimization in Finance
Description
Functions, data and other R code from the book ‘Numerical Methods and Optimization in Finance’. Comments/corrections/remarks/suggestions are very welcome (please contact the maintainer directly).
Details
The package contains implementations of several optimisation
heuristics: Differential Evolution (DEopt
), Genetic
Algorithms (GAopt
), (Stochastic) Local Search
(LSopt
), Particle Swarm (PSopt
),
Simuleated Annealing (SAopt
) and
Threshold Accepting (TAopt
). The term heuristic is meant
in the sense of generalpurpose optimisation method.
Dependencies: The package is completely written in R. A
number of packages are suggested, but they are not
strictly required when using the NMOF package, and
most of the package's functionality is available without
them. Specifically, package MASS is needed to run
the complete example for PSopt
and also in
one of the vignettes (PSlms
). Package
parallel is optional for functions
bracketing
, GAopt
,
gridSearch
and restartOpt
, and
may become an option for other functions. Package
quadprog is needed for a vignette
(TAportfolio
), some tests, and it may be used for
computing meanvariance efficient portfolios.
Package Rglpk is needed for function minCVaR
.
Package readxl is needed to process the raw data in function
Shiller
; package datetimeutils is
used by French
and Shiller
.
PMwR would be needed to run the examples
of the backtesting examples in the NMOF book.
Finally, packages RUnit and tinytest are needed to
run the tests in subdirectory ‘unitTests
’.
Version numbering: package versions are numbered in the
form majorminorpatch
. The patch level is
incremented with any published change in a version.
Minor version numbers are incremented when a
feature is added or an existing feature is substantially
revised. (Such changes will be reported in the NEWS file.)
The major version number will only be increased if
there were a new edition of the book.
The source code of the NMOF package is also hosted at https://github.com/enricoschumann/NMOF/. Updates to the package and new features are described at http://enricoschumann.net/notes/NMOF/.
Optimisation
There are functions for
Differential Evolution (DEopt
),
Genetic Algorithms (GAopt
),
(Stochastic) Local Search (LSopt
),
Simuleated Annealing (SAopt
),
Particle Swarm (SAopt
),
and Threshold Accepting (TAopt
).
The function restartOpt
helps with
running restarts of these methods;
also available are functions for
grid search (gridSearch
) and
greedy search (greedySearch
).
Pricing Financial Instruments
For options: See vanillaOptionEuropean
,
vanillaOptionAmerican
, putCallParity
.
For pricing methods that use the characteristic function, see
callCF
.
For bonds and bond futures: See vanillaBond
,
bundFuture
and xtContractValue
.
Simulation
Data
See bundData
, fundData
and
optionData
.
Author(s)
Enrico Schumann
Maintainer: Enrico Schumann <es@enricoschumann.net>
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
## Not run:
library("NMOF")
## overview
packageDescription("NMOF")
help(package = "NMOF")
## code from book
showExample("equations.R", edition = 1)
showExample("Heur")
## show NEWS file
news(Version >= "2.00", package = "NMOF")
## vignettes
vignette(package = "NMOF")
nss < vignette("DEnss", package = "NMOF")
print(nss)
edit(nss)
## _book_ websites
browseURL("http://nmof.net")
browseURL("http://enricoschumann.net/NMOF/")
## _package_ websites
browseURL("http://enricoschumann.net/R/packages/NMOF/")
browseURL("https://cran.rproject.org/package=NMOF")
browseURL("https://git.sr.ht/~enricoschumann/NMOF")
browseURL("https://github.com/enricoschumann/NMOF")
## unit tests
file.show(system.file("unitTests/test_results.txt", package = "NMOF"))
## End(Not run)
test.rep < readLines(system.file("unitTests/test_results.txt",
package = "NMOF"))
nt < gsub(".*\\(([09]+) checks?\\).*", "\\1",
test.rep[grep("\\(\\d+ checks?\\)", test.rep)])
message("Number of unit tests: ", sum(as.numeric(nt)))
ZeroBracketing
Description
Bracket the zeros (roots) of a univariate function
Usage
bracketing(fun, interval, ...,
lower = min(interval), upper = max(interval),
n = 20L,
method = c("loop", "vectorised", "multicore", "snow"),
mc.control = list(), cl = NULL)
Arguments
fun 
a univariate function; it will be called as 
interval 
a numeric vector, containing the endpoints of the interval to be searched 
... 
further arguments passed to 
lower 
lower endpoint. Ignored if 
upper 
upper endpoint. Ignored if 
n 
the number of function evaluations. Must be at least 2 (in which
case 
method 
can be 
mc.control 
a list containing settings that will be passed to 
cl 
default is 
Details
bracketing
evaluates fun
at equalspaced values of x
between (and including) lower
and upper
. If the sign of
fun
changes between two consecutive x
values,
bracketing
reports these two x
values as containing (‘bracketing’)
a root. There is no guarantee that there is only one root
within a reported interval. bracketing
will not narrow the chosen intervals.
The argument method
determines how fun
is
evaluated. Default is loop
. If method
is
"vectorised"
, fun
must be written such that it can be
evaluated for a vector x
(see Examples). If method
is
multicore
, function mclapply
from package parallel
is used. Further settings for mclapply
can be passed through
the list mc.control
. If multicore
is chosen but the
functionality is not available (eg, currently on Windows), then
method
will be set to loop
and a warning is issued. If
method
is snow
, function clusterApply
from
package parallel is used. In this case, the argument cl
must either be a cluster object (see the documentation of
clusterApply
) or an integer. If an integer, a cluster will be
set up via makeCluster(c(rep("localhost", cl)), type = "SOCK")
,
and stopCluster
is called when the function is exited. If
snow
is chosen but the package is not available or cl
is
not specified, then method
will be set to loop
and a
warning is issued. In case that cl
is a cluster object,
stopCluster
will not be called automatically.
Value
A numeric matrix with two columns, named lower and upper. Each row contains one interval that contains at least one root. If no roots were found, the matrix has zero rows.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
uniroot
(in package stats)
Examples
## Gilli/Maringer/Schumann (2011), p. 290
testFun < function(x)
cos(1/x^2)
bracketing(testFun, interval = c(0.3, 0.9), n = 26L)
bracketing(testFun, interval = c(0.3, 0.9), n = 26L, method = "vectorised")
German Government Bond Data
Description
A sample of data on 44 German government bonds. Contains ISIN, coupon, maturity and dirty price as of 20100531.
Usage
bundData
Format
bundData
is a list with three components: cfList
,
tmList
and bM
.
cfList
is list of 44 numeric vectors (the cash
flows). tmList
is a list of 44 character vectors (the payment dates)
formatted as YYYYMMDD. bM
is a
numeric vector with 44 elements (the dirty prices of the bonds).
Details
All prices are as of 31 May 2010. See chapter 14 in Gilli et al. (2011).
Source
The data was obtained from https://www.deutschefinanzagentur.de/en/ . The data is also freely available from the website of the Bundesbank https://www.bundesbank.de/en/ .
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
bundData
str(bundData)
## get ISINs of bonds
names(bundData$cfList)
## get a specific bond
thisBond < "DE0001135358"
data.frame(dates = as.Date(bundData$tmList[[thisBond]]),
payments = bundData$cfList[[thisBond]])
Theoretical Valuation of Euro Bund Future
Description
Compute theoretical prices of bund future.
Usage
bundFuture(clean, coupon, trade.date,
expiry.date, last.coupon.date,
r, cf)
bundFutureImpliedRate(future, clean, coupon,
trade.date, expiry.date,
last.coupon.date, cf)
Arguments
clean 
numeric: clean prices of CTD 
future 
numeric: price of future 
coupon 
numeric 
trade.date 

expiry.date 

last.coupon.date 

r 
numeric: 0.01 
cf 
numeric: conversion factor of CTD 
Details
bundFuture
computes the theoretical prices of the Bund Future,
given the prices of the cheapesttodeliver eligible government bond.
bundFutureImpliedRate
computes the implied refinancing rate.
Value
numeric
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
## BundFuture with expiry Sep 2017
## CTD: DE0001102408  0%, 15 Aug 2026
##
## On 21 August 2017, the CTD traded (clean) at 97.769
## the FGBL Sep 2017 closed at 164.44.
bundFuture(clean = 97.769, ## DE0001102408
coupon = 0,
trade.date = "2017821",
expiry.date = "20170907", ## Bund expiry
last.coupon.date = "20170815", ## last co
r = 0.0037,
cf = 0.594455) ## conversion factor (from Eurex website)
bundFutureImpliedRate(future = 164.44,
clean = 97.769,
coupon = 0,
trade.date = "2017821",
expiry.date = "20170907",
last.coupon.date = "20170815",
cf = 0.594455)
Price a PlainVanilla Call with the Characteristic Function
Description
Price a European plainvanilla call with the characteric function.
Usage
callCF(cf, S, X, tau, r, q = 0, ...,
implVol = FALSE, uniroot.control = list(), uniroot.info = FALSE)
cfBSM(om, S, tau, r, q, v)
cfMerton(om, S, tau, r, q, v, lambda, muJ, vJ)
cfBates(om, S, tau, r, q, v0, vT, rho, k, sigma, lambda, muJ, vJ)
cfHeston(om, S, tau, r, q, v0, vT, rho, k, sigma)
cfVG(om, S, tau, r, q, nu, theta, sigma)
Arguments
cf 
characteristic function 
S 
spot 
X 
strike 
tau 
time to maturity 
r 
the interest rate 
q 
the dividend rate 
... 
arguments passed to the characteristic function 
implVol 
logical: compute implied vol? 
uniroot.control 
A list. If there are elements named

uniroot.info 
logical; default is 
om 
a (usually complex) argument 
v0 
a numeric vector of length one 
vT 
a numeric vector of length one 
v 
a numeric vector of length one 
rho 
a numeric vector of length one 
k 
a numeric vector of length one 
sigma 
a numeric vector of length one 
lambda 
a numeric vector of length one 
muJ 
a numeric vector of length one 
vJ 
a numeric vector of length one 
nu 
a numeric vector of length one 
theta 
a numeric vector of length one 
Details
The function computes the value of a plain vanilla European call under
different models, using the representation of Bakshi/Madan. Put values
can be computed through put–call parity (see
putCallParity
).
If implVol
is TRUE
, the function will compute the
implied volatility necessary to obtain the same value under
Black–Scholes–Merton. The implied volatility is computed with
uniroot
from the stats package. The default search
interval is c(0.00001, 2)
; it can be changed through
uniroot.control
.
The function uses variances as inputs (not volatilities).
The function is not vectorised (but see the NMOF Manual for examples of how to efficiently price more than one option at once).
Value
Returns the value of the call (numeric) under the respective model or,
if implVol
is TRUE
, a list of the value and the implied
volatility. (If, in addition, uniroot.info
is TRUE
, the
information provided by uniroot
is also returned.)
Note
If implVol
is TRUE
, the function will return a list with
elements named value
and impliedVol
. Prior to version
0.263, the first element was named callPrice
.
Author(s)
Enrico Schumann
References
Bates, David S. (1996) Jumps and Stochastic Volatility: Exchange Rate Processes Implicit in Deutsche Mark Options. Review of Financial Studies 9 (1), 69–107.
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Heston, S.L. (1993) A ClosedForm Solution for Options with Stochastic Volatility with Applications to Bonds and Currency options. Review of Financial Studies 6 (2), 327–343.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
S < 100; X < 100; tau < 1
r < 0.02; q < 0.08
v0 < 0.2^2 ## variance, not volatility
vT < 0.2^2 ## variance, not volatility
v < vT
rho < 0.3; k < .2
sigma < 0.3
## jump parameters (Merton and Bates)
lambda < 0.1
muJ < 0.2
vJ < 0.1^2
## get Heston price and BSM implied volatility
callHestoncf(S, X, tau, r, q, v0, vT, rho, k, sigma, implVol = FALSE)
callCF(cf = cfHeston, S=S, X=X, tau=tau, r=r, q = q,
v0 = v0, vT = vT, rho = rho, k = k, sigma = sigma, implVol = FALSE)
## BlackScholesMerton
callCF(cf = cfBSM, S=S, X=X, tau = tau, r = r, q = q,
v = v, implVol = TRUE)
## Bates
callCF(cf = cfBates, S = S, X = X, tau = tau, r = r, q = q,
v0 = v0, vT = vT, rho = rho, k = k, sigma = sigma,
lambda = lambda, muJ = muJ, vJ = vJ, implVol = FALSE)
## Merton
callCF(cf = cfMerton, S = S, X = X, tau = tau, r = r, q = q,
v = v, lambda = lambda, muJ = muJ, vJ = vJ, implVol = FALSE)
## variance gamma
nu < 0.1; theta < 0.1; sigma < 0.15
callCF(cf = cfVG, S = S, X = X, tau = tau, r = r, q = q,
nu = nu, theta = theta, sigma = sigma, implVol = FALSE)
Price of a European Call under the Heston Model
Description
Computes the price of a European Call under the Heston model (and the equivalent Black–Scholes–Merton volatility)
Usage
callHestoncf(S, X, tau, r, q, v0, vT, rho, k, sigma, implVol = FALSE,
...,
uniroot.control = list(), uniroot.info = FALSE)
Arguments
S 
current stock price 
X 
strike price 
tau 
time to maturity 
r 
riskfree rate 
q 
dividend rate 
v0 
current variance 
vT 
longrun variance (theta in Heston's paper) 
rho 
correlation between spot and variance 
k 
speed of meanreversion (kappa in Heston's paper) 
sigma 
volatility of variance. A value smaller than 0.01 is replaced with 0.01. 
implVol 
compute equivalent Black–Scholes–Merton
volatility? Default is 
... 
named arguments, passed to 
uniroot.control 
A list. If there are elements named

uniroot.info 
logical; default is 
Details
The function computes the value of a plain vanilla European call under the Heston model. Put values can be computed through put–callparity.
If implVol
is TRUE
, the function will
compute the implied volatility necessary to obtain the
same price under Black–Scholes–Merton. The implied
volatility is computed with uniroot
from
the stats package (the default search interval is
c(0.00001, 2)
; it can be changed through
uniroot.control
).
Note that the function takes variances as inputs (not volatilities).
Value
Returns the value of the call (numeric) under the Heston
model or, if implVol
is TRUE
, a list of the
value and the implied volatility. If uniroot.info
is TRUE
, then instead of only the computed
volatility, the complete output of uniroot
is included in the result.
Note
If implVol
is TRUE
, the function will
return a list with elements named value
and
impliedVol
. Prior to version 0.263, the first
element was named callPrice
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Heston, S.L. (1993) A ClosedForm Solution for Options with Stochastic Volatility with Applications to Bonds and Currency options. Review of Financial Studies 6(2), 327–343.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
S < 100; X < 100; tau < 1; r < 0.02; q < 0.01
v0 < 0.2^2 ## variance, not volatility
vT < 0.2^2 ## variance, not volatility
rho < 0.7; k < 0.2; sigma < 0.5
## get Heston price and BSM implied volatility
result < callHestoncf(S = S, X = X, tau = tau, r = r, q = q,
v0 = v0, vT = vT, rho = rho, k = k,
sigma = sigma, implVol = TRUE)
## Heston price
result[[1L]]
## price BSM with implied volatility
vol < result[[2L]]
d1 < (log(S/X) + (r  q + vol^2 / 2)*tau) / (vol*sqrt(tau))
d2 < d1  vol*sqrt(tau)
callBSM < S * exp(q * tau) * pnorm(d1) 
X * exp(r * tau) * pnorm(d2)
callBSM ## should be (about) the same as result[[1L]]
Price of a European Call under Merton's Jump–Diffusion Model
Description
Computes the price of a European Call under Merton's jump–diffusion model (and the equivalent Black–Scholes–Merton volatility)
Usage
callMerton(S, X, tau, r, q, v, lambda, muJ, vJ, N, implVol = FALSE)
Arguments
S 
current stock price 
X 
strike price 
tau 
time to maturity 
r 
riskfree rate 
q 
dividend rate 
v 
variance 
lambda 
jump intensity 
muJ 
mean jumpsize 
vJ 
variance of log jumpsize 
N 
The number of jumps. See Details. 
implVol 
compute equivalent Black–Scholes–Merton volatility?
Default is 
Details
The function computes the value of a plainvanilla European call under
Merton's jump–diffusion model. Put values can be computed through
put–callparity (see putCallParity
). If implVol
is TRUE
, the function also computes the implied volatility
necessary to obtain the same price under Black–Scholes–Merton. The
implied volatility is computed with uniroot
from the
stats package.
Note that the function takes variances as inputs (not volatilities).
The number of jumps N
typically can be set 10 or 20. (Just try to
increase N
and see how the results change.)
Value
Returns the value of the call (numeric) or, if implVol
is
TRUE
, a list of the value and the implied volatility.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Merton, R.C. (1976) Option Pricing when Underlying Stock Returns are Discontinuous. Journal of Financial Economics 3(1–2), 125–144.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
S < 100; X < 100; tau < 1
r < 0.0075; q < 0.00
v < 0.2^2
lambda < 1; muJ < 0.2; vJ < 0.6^2
N < 20
## jumps can make a difference
callMerton(S, X, tau, r, q, v, lambda, muJ, vJ, N, implVol = TRUE)
callCF(cf = cfMerton, S = S, X = X, tau = tau, r = r, q = q,
v = v, lambda = lambda, muJ = muJ, vJ = vJ, implVol = TRUE)
vanillaOptionEuropean(S,X,tau,r,q,v, greeks = FALSE)
lambda < 0 ## no jumps
callMerton(S, X, tau, r, q, v, lambda, muJ, vJ, N, implVol = FALSE)
vanillaOptionEuropean(S,X,tau,r,q,v, greeks = FALSE)
lambda < 1; muJ < 0; vJ < 0.0^2 ## no jumps, either
callMerton(S, X, tau, r, q, v, lambda, muJ, vJ, N, implVol = FALSE)
vanillaOptionEuropean(S,X,tau,r,q,v, greeks = FALSE)
Fullrank Column Subset
Description
Select a fullrank subset of columns of a matrix.
Usage
colSubset(x)
Arguments
x 
a numeric matrix 
Details
Uses qr
.
Value
A list:
columns 
indices of columns 
multiplier 
a matrix 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
nc < 3 ## columns
nr < 10 ## rows
M < array(rnorm(nr * nc), dim = c(nr, nc))
C < array(0.5, dim = c(nc, nc))
diag(C) < 1
M < M %*% chol(C)
M < M[ ,c(1,1,1,2,3)]
M
(tmp < colSubset(M))
C < cor(M[ ,tmp$columns])
nc < ncol(C)
nr < 100
X < array(rnorm(nr*nc), dim = c(nr, nc))
X < X %*% chol(C)
X < X %*% tmp$multiplier
head(X)
cor(X)
ConstantProportion Portfolio Insurance
Description
Simulate constantproportion portfolio insurance (CPPI) for a given price path.
Usage
CPPI(S, multiplier, floor, r, tau = 1, gap = 1)
Arguments
S 
numeric: price path of risky asset 
multiplier 
numeric 
floor 
numeric: a percentage, should be smaller than 1 
r 
numeric: interest rate (per time period tau) 
tau 
numeric: time periods 
gap 
numeric: how often to rebalance. 1 means every timestep, 2 means every second timestep, and so on. 
Details
Based on Dietmar Maringer's MATLAB code (function CPPIgap, Listing 9.1).
See Gilli, Maringer and Schumann, 2011, chapter 9.
Value
A list:
V 
normalised value (always starts at 1) 
C 
cushion 
B 
bond investment 
F 
floor 
E 
exposure 
N 
units of risky asset 
S 
price path 
Author(s)
Original MATLAB code: Dietmar Maringer. R implementation: Enrico Schumann.
References
Chapter 9 of Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
tau < 2
S < gbm(npaths = 1, timesteps = tau*256,
r = 0.02, v = 0.2^2, tau = tau, S0 = 100)
## rebalancing every day
sol < CPPI(S, multiplier = 5, floor = 0.9, r = 0.01,
tau = tau, gap = 1)
par(mfrow = c(3,1), mar = c(3,3,1,1))
plot(0:(length(S)1), S, type = "s", main = "stock price")
plot(0:(length(S)1), sol$V, type = "s", main = "value")
plot(0:(length(S)1), 100*sol$E/sol$V, type = "s",
main = "% invested in risky asset")
## rebalancing every 5th day
sol < CPPI(S, multiplier = 5, floor = 0.9, r = 0.01,
tau = tau, gap = 5)
par(mfrow = c(3,1), mar = c(3,3,1,1))
plot(0:(length(S)1), S, type = "s", main = "stock price")
plot(0:(length(S)1), sol$V, type = "s", main = "value")
plot(0:(length(S)1), 100*sol$E/sol$V, type = "s",
main = "% invested in risky asset")
Optimisation with Differential Evolution
Description
The function implements the standard Differential Evolution algorithm.
Usage
DEopt(OF, algo = list(), ...)
Arguments
OF 
The objective function, to be minimised. See Details. 
algo 
A list with the settings for algorithm. See Details and Examples. 
... 
Other pieces of data required to evaluate the objective function. See Details and Examples. 
Details
The function implements the standard Differential Evolution (no jittering or other features). Differential Evolution (DE) is a populationbased optimisation heuristic proposed by Storn and Price (1997). DE evolves several solutions (collected in the ‘population’) over a number of iterations (‘generations’). In a given generation, new solutions are created and evaluated; better solutions replace inferior ones in the population. Finally, the best solution of the population is returned. See the references for more details on the mechanisms.
To allow for constraints, the evaluation works as follows: after a new
solution is created, it is (i) repaired, (ii) evaluated through the
objective function, (iii) penalised. Step (ii) is done by a call to
OF
; steps (i) and (iii) by calls to algo$repair
and
algo$pen
. Step (i) and (iii) are optional, so the respective
functions default to NULL
. A penalty is a positive number added
to the ‘clean’ objective function value, so it can also be
directly written in the OF
. Writing a separate penalty function
is often clearer; it can be more efficient if either only the objective
function or only the penalty function can be vectorised. (Constraints
can also be added without these mechanisms. Solutions that violate
constraints can, for instance, be mapped to feasible solutions, but
without actually changing them. See Maringer and Oyewumi, 2007, for an
example.)
Conceptually, DE consists of two loops: one loop across the
generations and, in any given generation, one loop across the solutions.
DEopt
indeed uses, as the default, two loops. But it does not
matter in what order the solutions are evaluated (or repaired or
penalised), so the second loop can be vectorised. This is controlled by
the variables algo$loopOF
, algo$loopRepair
and
algo$loopPen
, which all default to TRUE
. Examples are
given in the vignettes and in the book. The respective
algo$loopFun
must then be set to FALSE
.
All objects that are passed through ...
will be passed to the
objective function, to the repair function and to the penalty function.
The list algo
collects the the settings for the
algorithm. Strictly necessary are only min
and max
(to
initialise the population). Here are all possible arguments:
CR
probability for crossover. Defaults to 0.9. Using default settings may not be a good idea.
F
The step size. Typically a numeric vector of length one; default is 0.5. Using default settings may not be a good idea. (
F
can also be a vector with different values for each decision variable.)nP
population size. Defaults to 50. Using default settings may not be a good idea.
nG
number of generations. Defaults to 300. Using default settings may not be a good idea.
min
,max
vectors of minimum and maximum parameter values. The vectors
min
andmax
are used to determine the dimension of the problem and to randomly initialise the population. Per default, they are no constraints: a solution may well be outside these limits. Only ifalgo$minmaxConstr
isTRUE
will the algorithm repair solutions outside themin
andmax
range.minmaxConstr
if
TRUE
,algo$min
andalgo$max
are considered constraints. Default isFALSE
.pen
a penalty function. Default is
NULL
(no penalty).initP
optional: the initial population. A matrix of size
length(algo$min)
timesalgo$nP
, or a function that creates such a matrix. If a function, it should take no arguments.repair
a repair function. Default is
NULL
(no repairing).loopOF
logical. Should the
OF
be evaluated through a loop? Defaults toTRUE
.loopPen
logical. Should the penalty function (if specified) be evaluated through a loop? Defaults to
TRUE
.loopRepair
logical. Should the repair function (if specified) be evaluated through a loop? Defaults to
TRUE
.printDetail
If
TRUE
(the default), information is printed. If an integeri
greater then one, information is printed at veryi
th generation.printBar
If
TRUE
(the default), atxtProgressBar
is printed.storeF
if
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
default is
FALSE
. IfTRUE
, the solutions (ie, decision variables) in every generation are stored and returned as a listP
in listxlist
(see Value section below). To check, for instance, the solutions at the end of thei
th generation, retrievexlist[[c(1L, i)]]
. This will be a matrix of sizelength(algo$min)
timesalgo$nP
. (To be consistent with other functions,xlist
is itself a list. In the case ofDEopt
, it contains just one element.)classify
Logical; default is
FALSE
. IfTRUE
, the result will have a class attributeTAopt
attached. This feature is experimental: the supported methods may change without warning.drop

If
FALSE
(the default), the dimension is not dropped from a single solution when it is passed to a function. (That is, the function will receive a singlecolumn matrix.)
Value
A list:
xbest 
the solution (the best member of the population), which is a numeric vector 
OFvalue 
objective function value of best solution 
popF 
a vector. The objective function values in the final population. 
Fmat 
if 
xlist 
if 
initial.state 
the value of 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Maringer, D. and Oyewumi, O. (2007). Index Tracking with Constrained Portfolios. Intelligent Systems in Accounting, Finance and Management, 15(1), pp. 57–71.
Schumann, E. (2012) Remarks on 'A comparison of some heuristic optimization methods'. http://enricoschumann.net/R/remarks.htm
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Storn, R., and Price, K. (1997) Differential Evolution – a Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. Journal of Global Optimization, 11(4), pp. 341–359.
See Also
Examples
## Example 1: Trefethen's 100digit challenge (problem 4)
## http://people.maths.ox.ac.uk/trefethen/hundred.html
OF < tfTrefethen ### see ?testFunctions
algo < list(nP = 50L, ### population size
nG = 300L, ### number of generations
F = 0.6, ### step size
CR = 0.9, ### prob of crossover
min = c(10, 10), ### range for initial population
max = c( 10, 10))
sol < DEopt(OF = OF, algo = algo)
## correct answer: 3.30686864747523
format(sol$OFvalue, digits = 12)
## check convergence of population
sd(sol$popF)
ts.plot(sol$Fmat, xlab = "generations", ylab = "OF")
## Example 2: vectorising the evaluation of the population
OF < tfRosenbrock ### see ?testFunctions
size < 3L ### define dimension
x < rep.int(1, size) ### the known solution ...
OF(x) ### ... should give zero
algo < list(printBar = FALSE,
nP = 30L,
nG = 300L,
F = 0.6,
CR = 0.9,
min = rep(100, size),
max = rep( 100, size))
## run DEopt
(t1 < system.time(sol < DEopt(OF = OF, algo = algo)))
sol$xbest
sol$OFvalue ### should be zero (with luck)
## a vectorised Rosenbrock function: works only with a *matrix* x
OF2 < function(x) {
n < dim(x)[1L]
xi < x[seq_len(n  1L), ]
colSums(100 * (x[2L:n, ]  xi * xi)^2 + (1  xi)^2)
}
## random solutions (every column of 'x' is one solution)
x < matrix(rnorm(size * algo$nP), size, algo$nP)
all.equal(OF2(x)[1:3],
c(OF(x[ ,1L]), OF(x[ ,2L]), OF(x[ ,3L])))
## run DEopt and compare computing time
algo$loopOF < FALSE
(t2 < system.time(sol2 < DEopt(OF = OF2, algo = algo)))
sol2$xbest
sol2$OFvalue ### should be zero (with luck)
t1[[3L]]/t2[[3L]] ### speedup
Diversification Ratio
Description
Compute the diversification ratio of a portfolio.
Usage
divRatio(w, var)
Arguments
w 
numeric: a vector of weights 
var 
numeric matrix: the variance–covariance matrix 
Details
The function provides an efficient implementation of the diversification ratio, suitable for optimisation.
Value
a numeric vector of length one
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Yves Choueifaty and Yves Coignard (2008) Toward Maximum Diversification. Journal of Portfolio Management 35(1), 40–51.
See Also
pm
, drawdown
Examples
na < 10 ## number of assets
rho < 0.5 ## correlation
v_min < 0.2 ## minimum vol
v_max < 0.4 ## maximum vol
## set up a covariance matrix S
C < array(rho, dim = c(na,na))
diag(C) < 1
vols < seq(v_min, v_max, length.out = na)
S < outer(vols, vols) * C
w < rep(1/na, na) ## weights
divRatio(w, S)
Drawdown
Description
Compute the drawdown of a time series.
Usage
drawdown(v, relative = TRUE, summary = TRUE)
Arguments
v 
a price series (a numeric vector) 
relative 
if 
summary 
if 
Details
The drawdown at position t of a time series v is the difference between the highest peak that was reached before t and the current value. If the current value represents a new high, the drawdown is zero.
Value
If summary
is FALSE
, a vector of the same length as
v
. If summary
is TRUE
, a list
maximum 
maximum drawdown 
high 
the max of 
high.position 
position of 
low 
the min of 
low.position 
position of 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
v < cumprod(1 + rnorm(20) * 0.02)
drawdown(v)
Computing Prices of European Calls with a Binomial Tree
Description
Computes the fair value of a European Call with the binomial tree of Cox, Ross and Rubinstein.
Usage
EuropeanCall(S0, X, r, tau, sigma, M = 101)
EuropeanCallBE(S0, X, r, tau, sigma, M = 101)
Arguments
S0 
current stock price 
X 
strike price 
r 
riskfree rate 
tau 
time to maturity 
sigma 
volatility 
M 
number of time steps 
Details
Prices a European Call with the tree approach of Cox, Ross, Rubinstein.
The algorithm in EuropeanCallBE
does not construct and traverse a
tree, but computes the terminal prices via a binomial expansion (see
Higham, 2002, and Chapter 5 in Gilli/Maringer/Schumann, 2011).
Value
Returns the value of the call (numeric
).
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
M. Gilli and Schumann, E. (2009) Implementing Binomial Trees. COMISEF Working Paper Series No. 008. http://enricoschumann.net/COMISEF/wps008.pdf
Higham, D. (2002) Nine Ways to Implement the Binomial Method for Option Valuation in MATLAB. SIAM Review, 44(4), pp. 661–677. doi:10.1137/S0036144501393266 .
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## price
EuropeanCall( S0 = 100, X = 100, r = 0.02, tau = 1, sigma = 0.20, M = 50)
EuropeanCallBE(S0 = 100, X = 100, r = 0.02, tau = 1, sigma = 0.20, M = 50)
## a Greek: delta
h < 1e8
C1 < EuropeanCall(S0 = 100 + h, X = 100, r = 0.02, tau = 1,
sigma = 0.20, M = 50)
C2 < EuropeanCall(S0 = 100 , X = 100, r = 0.02, tau = 1,
sigma = 0.20, M = 50)
(C1  C2) / h
Download Datasets from Kenneth French's Data Library
Description
Download datasets from Kenneth French's Data Library.
Usage
French(dest.dir,
dataset = "FF_Research_Data_Factors_CSV.zip",
weighting = "value", frequency = "monthly",
price.series = FALSE, na.rm = FALSE,
adjust.frequency = TRUE)
Arguments
dest.dir 
character: a path to a directory 
dataset 
a character string: the CSV file name. Also
supported are the keywords ‘ 
weighting 
a character string: 
frequency 
a character string: 
price.series 
logical: convert the returns series into prices series? 
na.rm 
logical: remove missing values in the calculation of price series? 
adjust.frequency 
logical: if 
Details
The function downloads data provided by Kenneth
French at
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.
The download file gets a date prefix (current date in
format YYYYMMDD
) and is stored in directory
dest.dir
. Before any download is attempted,
the function checks whether a file with today's
prefix exist in dest.dir
; if yes, the file is
used.
In the original data files, missing values are
coded as 99
or similar. These
numeric values are replaced by NA
.
Calling the function without any arguments will print the names of the supported datasets (and return them insivibly).
Value
A data.frame
, with contents depending on the
particular dataset. If the download failes, the function
evaluates to NULL
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## list all supported files
French()
## fetch names of files from Kenneth French's website
try({
txt < readLines(paste0("https://mba.tuck.dartmouth.edu/pages/",
"faculty/ken.french/data_library.html"))
csv < txt[grep("ftp/.*CSV.zip", txt, ignore.case = TRUE)]
gsub(".*ftp/(.*?CSV.zip).*", "\1", csv, ignore.case = TRUE)
})
## Not run:
archive.dir < "~/Downloads/French"
if (!dir.exists(archive.dir))
dir.create(archive.dir)
French(archive.dir, "FF_Research_Data_Factors_CSV.zip")
## End(Not run)
Mutual Fund Returns
Description
A matrix of 500 rows (return scenarios) and 200 columns (mutual funds). The elements in the matrix are weekly returns.
Usage
fundData
Format
A plain numeric matrix.
Details
The scenarios were created with a bootstrapping technique. The data set is only meant to provide example data on which to test algorithms.
Source
Schumann, E. (2010) Essays on Practical Financial Optimisation, (chapter 4), PhD thesis, University of Geneva.
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
apply(fundData, 2, summary)
Optimisation with a Genetic Algorithm
Description
A simple Genetic Algorithm for minimising a function.
Usage
GAopt (OF, algo = list(), ...)
Arguments
OF 
The objective function, to be minimised. See Details. 
algo 
A list with the settings for algorithm. See Details and Examples. 
... 
Other pieces of data required to evaluate the objective function. See Details and Examples. 
Details
The function implements a simple Genetic Algorithm (GA). A
GA evolves a collection of solutions (the socalled
population), all of which are coded as vectors containing only zeros
and ones. (In GAopt
, solutions are of mode logical
.)
The algorithm starts with randomlychosen or usersupplied population
and aims to iteratively improve this population by mixing solutions
and by switching single bits in solutions, both at random. In each
iteration, such randomlychanged solutions are compared with the
original population and better solutions replace inferior
ones. In GAopt
, the population size is kept constant.
GA language: iterations are called generations; new solutions
are called offspring or children (and the existing solutions, from which
the children are created, are parents); the objective function is called
a fitness function; mixing solutions is a crossover; and randomly
changing solutions is called mutation. The choice which solutions remain in
the population and which ones are discarded is called selection. In
GAopt
, selection is pairwise: a given child is compared with a
given parent; the better of the two is kept. In this way, the best
solution is automatically retained in the population.
To allow for constraints, the evaluation works as follows: after new
solutions are created, they are (i) repaired, (ii) evaluated through the
objective function, (iii) penalised. Step (ii) is done by a call to
OF
; steps (i) and (iii) by calls to algo$repair
and
algo$pen
. Step (i) and (iii) are optional, so the respective
functions default to NULL
. A penalty can also be directly written
in the OF
, since it amounts to a positive number added to the
‘clean’ objective function value; but a separate function is
often clearer. A separate penalty function is advantagous if either only
the objective function or only the penalty function can be vectorised.
Conceptually a GA consists of two loops: one loop across the
generations and, in any given generation, one loop across the solutions.
This is the default, controlled by the variables algo$loopOF
,
algo$loopRepair
and algo$loopPen
, which all default to
TRUE
. But it does not matter in what order the solutions are
evaluated (or repaired or penalised), so the second loop can be
vectorised. The respective algo$loopFun
must then be set to
FALSE
. (See also the examples for DEopt
and
PSopt
.)
The evaluation of the objective function in a given generation can even
be distributed. For this, an argument algo$methodOF
needs to be
set; see below for details (and Schumann, 2011, for examples).
All objects that are passed through ...
will be passed to the
objective function, to the repair function and to the penalty function.
The list algo
contains the following items:
nB
number of bits per solution. Must be specified.
nP
population size. Defaults to 50. Using default settings may not be a good idea.
nG
number of iterations (‘generations’). Defaults to 300. Using default settings may not be a good idea.
crossover
The crossover method. Default is
"onePoint"
; also possible is “uniform”.prob
The probability for switching a single bit. Defaults to 0.01; typically a small number.
pen
a penalty function. Default is
NULL
(no penalty).repair
a repair function. Default is
NULL
(no repairing).initP
optional: the initial population. A logical matrix of size
length(algo$nB)
timesalgo$nP
, or a function that creates such a matrix. If a function, it must take no arguments. Ifmode(mP)
is notlogical
, thenstorage.mode(mP)
will be tried (and a warning will be issued).loopOF
logical. Should the
OF
be evaluated through a loop? Defaults toTRUE
.loopPen
logical. Should the penalty function (if specified) be evaluated through a loop? Defaults to
TRUE
.loopRepair
logical. Should the repair function (if specified) be evaluated through a loop? Defaults to
TRUE
.methodOF
loop
(the default),vectorised
,snow
ormulticore
. Settingvectorised
is equivalent to havingalgo$loopOF
set toFALSE
(andmethodOF
overridesloopOF
).snow
andmulticore
use functionsclusterApply
andmclapply
, respectively. Forsnow
, an objectalgo$cl
needs to be specified (see below). Formulticore
, optional arguments can be passed throughalgo$mc.control
(see below).cl
a cluster object or the number of cores. See documentation of package
parallel
.mc.control
a list of named elements; optional settings for
mclapply
(for instance,list(mc.set.seed = FALSE)
)printDetail
If
TRUE
(the default), information is printed.printBar
If
TRUE
(the default), atxtProgressBar
is printed.storeF
If
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
If
TRUE
, the solutions (ie, binary strings) in every generation are stored and returned as a listP
in listxlist
(see Value section below). To check, for instance, the solutions at the end of thei
th generation, retrievexlist[[c(1L, i)]]
. This will be a matrix of sizealgo$nB
timesalgo$nP
.classify
Logical; default is
FALSE
. IfTRUE
, the result will have a class attributeTAopt
attached. This feature is experimental: the supported methods may change without warning.
Value
A list:
xbest 
the solution (the best member of the population) 
OFvalue 
objective function value of best solution 
popF 
a vector. The objective function values in the final population. 
Fmat 
if 
xlist 
if 
initial.state 
the value of 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## a *very* simple problem (why?):
## match a binary (logical) string y
size < 20L ### the length of the string
OF < function(x, y) sum(x != y)
y < runif(size) > 0.5
x < runif(size) > 0.5
OF(y, y) ### the optimum value is zero
OF(x, y)
algo < list(nB = size, nP = 20L, nG = 100L, prob = 0.002)
sol < GAopt(OF, algo = algo, y = y)
## show differences (if any: marked by a '^')
cat(as.integer(y), "\n", as.integer(sol$xbest), "\n",
ifelse(y == sol$xbest , " ", "^"), "\n", sep = "")
algo$nP < 3L ### that shouldn't work so well
sol2 < GAopt(OF, algo = algo, y = y)
## show differences (if any: marked by a '^')
cat(as.integer(y), "\n", as.integer(sol2$xbest), "\n",
ifelse(y == sol2$xbest , " ", "^"), "\n", sep = "")
Greedy Search
Description
Greedy Search
Usage
greedySearch(OF, algo, ...)
Arguments
OF 
The objective function, to be minimised. Its first
argument needs to be a solution; 
algo 
List of settings. See Details. 
... 
Other variables to be passed to the objective function and to the neighbourhood function. See Details. 
Details
A greedy search works starts at a provided initial solution (called the current solution) and searches a defined neighbourhood for the best possible solution. If this best neighbour is not better than the current solution, the search stops. Otherwise, the best neighbour becomes the current solution, and the search is repeated.
Value
A list:
xbest 
best solution found. 
OFvalue 
objective function value associated with best solution. 
Fmat 
a matrix with two
columns. 
xlist 
a list 
initial.state 
the value of

x0 
the initial solution 
iterations 
the number of iterations after which the search stopped 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
LSopt
Examples
na < 100
inc < 5
R < randomReturns(na = na,
ns = 1000,
sd = seq(0.01, 0.02, length.out = 100),
rho = 0.5)
S < cov(R)
OF < function(x, S, ...) {
w < 1/sum(x)
sum(w * w * S[x, x])
}
x < logical(na)
x[1:inc] < TRUE
all.neighbours < function(x, ...) {
true < which( x)
false < which(!x)
ans < list()
for (i in true) {
for (j in false) {
ans1 < x
ans1[i] < !x[i]
ans1[j] < !x[j]
ans < c(ans, list(ans1))
}
}
ans
}
algo < list(loopOF = TRUE,
maxit = 1000,
all.neighbours = all.neighbours,
x0 = x)
system.time(sol.gs < greedySearch(OF, algo = algo, S = S))
sqrt(sol.gs$OFvalue)
Grid Search
Description
Evaluate a function for a given list of arguments.
Usage
gridSearch(fun, levels, ..., lower, upper, npar = 1L, n = 5L,
printDetail = TRUE,
method = NULL,
mc.control = list(), cl = NULL,
keepNames = FALSE, asList = FALSE)
Arguments
fun 
a function of the form 
levels 
a list of levels for the arguments (see Examples) 
... 
objects passed to 
lower 
a numeric vector. Ignored if levels are explicitly specified. 
upper 
a numeric vector. Ignored if levels are explicitly specified. 
npar 
the number of parameters. Must be supplied if 
n 
the number of levels. Default is 5. Ignored if levels are explicitly specified. 
printDetail 
print information on the number of objective function evaluations 
method 
can be 
mc.control 
a list containing settings that will be passed to 
cl 
default is 
keepNames 

asList 
does 
Details
A grid search can be used to find ‘good’ parameter values for a
function. In principle, a grid search has an obvious deficiency: as
the length of x
(the first argument to fun
) increases,
the number of necessary function evaluations grows exponentially. Note
that gridSearch
will not warn about an unreasonable number of
function evaluations, but if printDetail
is TRUE
it will
print the required number of function evaluations.
In practice, grid search is often better than its reputation. If a function takes only a few parameters, it is often a reasonable approach to find ‘good’ parameter values.
The function uses the mechanism of expand.grid
to create
the list of parameter combinations for which fun
is evaluated; it
calls lapply
to evaluate fun
if
method == "loop"
(the default).
If method
is multicore
, then function mclapply
from package parallel is used. Further settings for
mclapply
can be passed through the list mc.control
. If
multicore
is chosen but the functionality is not available,
then method
will be set to loop
and a warning is
issued. If method == "snow"
, the function clusterApply
from package parallel is used. In this case, the argument cl
must either be a cluster object (see the documentation of
clusterApply
) or an integer. If an integer, a cluster will be
set up via makeCluster(c(rep("localhost", cl)), type = "SOCK")
(and stopCluster
is called when the function is exited). If
snow
is chosen but not available or cl
is not specified,
then method
will be set to loop
and a warning is issued.
Value
A list.
minfun 
the minimum of 
minlevels 
the levels that give this minimum. 
values 
a list. All the function values of 
levels 
a list. All the levels for which 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
testFun < function(x)
x[1L] + x[2L]^2
sol < gridSearch(fun = testFun, levels = list(1:2, c(2, 3, 5)))
sol$minfun
sol$minlevels
## specify all levels
levels < list(a = 1:2, b = 1:3)
res < gridSearch(testFun, levels)
res$minfun
sol$minlevels
## specify lower, upper and npar
lower < 1; upper < 3; npar < 2
res < gridSearch(testFun, lower = lower, upper = upper, npar = npar)
res$minfun
sol$minlevels
## specify lower, upper, npar and n
lower < 1; upper < 3; npar < 2; n < 4
res < gridSearch(testFun, lower = lower, upper = upper, npar = npar, n = n)
res$minfun
sol$minlevels
## specify lower, upper and n
lower < c(1,1); upper < c(3,3); n < 4
res < gridSearch(testFun, lower = lower, upper = upper, n = n)
res$minfun
sol$minlevels
## specify lower, upper (autoexpanded) and n
lower < c(1,1); upper < 3; n < 4
res < gridSearch(testFun, lower = lower, upper = upper, n = n)
res$minfun
sol$minlevels
## nonnumeric inputs
test_fun < function(x) {
(length(x$S) + x$N1 + x$N2)
}
ans < gridSearch(test_fun,
levels = list(S = list("a", c("a", "b"), c("a", "b", "c")),
N1 = 1:5,
N2 = 101:105),
asList = TRUE, keepNames = TRUE)
ans$minlevels
## $S
## [1] "a" "b" "c"
##
## $N1
## [1] 5
##
## $N2
## [1] 105
LocalSearch Information
Description
The function can be called from the objective and neighbourhood
function during a run of LSopt
; it provides information
such as the current iteration.
Usage
LS.info(n = 0L)
Arguments
n 
generational offset; see Details. 
Details
This function is still experimental.
The function can be called in the neighbourhood function or the
objective function during a run of LSopt
. It evaluates
to a list with the state of the optimisation run, such as the current
iteration.
LS.info
relies on parent.frame
to retrieve its
information. If the function is called within another function in the
neighbourhood or objective function, the argument n
needs to be
increased.
Value
A list
iteration 
current iteration 
step 
same as ‘iteration’ 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## MINIMAL EXAMPLE for LSopt
## objective function evaluates to a constant
fun < function(x)
0
## neighbourhood function does not even change the solution,
## but it reports information
nb < function(x) {
tmp < LS.info()
cat("current iteration ", tmp$iteration, "\n")
x
}
## run LS
algo < list(nS = 5,
x0 = rep(0, 5),
neighbour = nb,
printBar = FALSE)
ignore < LSopt(fun, algo)
Stochastic Local Search
Description
Performs a simple stochastic Local Search.
Usage
LSopt(OF, algo = list(), ...)
Arguments
OF 
The objective function, to be minimised. Its first argument needs to
be a solution; 
algo 
List of settings. See Details. 
... 
Other variables to be passed to the objective function and to the neighbourhood function. See Details. 
Details
Local Search (LS) changes an initial solution for a number
of times, accepting only such changes that lead to an improvement in
solution quality (as measured by the objective function OF
).
More specifically, in each iteration, a current solution xc
is
changed through a function algo$neighbour
. This function takes
xc
as an argument and returns a new solution xn
. If
xn
is not worse than xc
, ie, if
OF(xn,...)<=OF(xc,...)
, then xn
replaces xc
.
The list algo
contains the following items:
nS
The number of steps. The default is 1000; but this setting depends very much on the problem.
nI

Total number of iterations, with default
NULL
. If specified, it will overridenS
. The option is provided to makes it easier to compare and switch between functionsLSopt
,TAopt
andSAopt
. x0
The initial solution. This can be a function; it will then be called once without arguments to compute an initial solution, ie,
x0 < algo$x0()
. This can be useful whenLSopt
is called in a loop of restarts and each restart is to have its own starting value.neighbour
The neighbourhood function, called as
neighbour(x, ...)
. Its first argument must be a solutionx
; it must return a changed solution.printDetail
If
TRUE
(the default), information is printed. If an integeri
greater then one, information is printed at veryi
th step.printBar
If
TRUE
(the default), atxtProgressBar
(from package utils) is printed). The progress bar is not shown ifprintDetail
is an integer greater than 1.storeF
if
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
default is
FALSE
. IfTRUE
, the solutions (ie, decision variables) in every generation are stored and returned in listxlist
(see Value section below). To check, for instance, the current solution at the end of thei
th generation, retrievexlist[[c(2L, i)]]
.OF.target
Numeric; when specified, the algorithm will stop when an objectivefunction value as low as
OF.target
(or lower) is achieved. This is useful when an optimal objectivefunction value is known: the algorithm will then stop and not waste time searching for a better solution.
At the minimum, algo
needs to contain an initial solution
x0
and a neighbour
function.
LS works on solutions through the functions neighbour
and OF
, which are specified by the user. Thus, a solution need
not be a numeric vector, but can be any other data structure as well
(eg, a list or a matrix).
To run silently (except for warnings and errors),
algo$printDetail
and algo$printBar
must be FALSE
.
Value
A list:
xbest 
best solution found. 
OFvalue 
objective function value associated with best solution. 
Fmat 
a matrix with two columns. 
xlist 
if 
initial.state 
the value of 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
TAopt
, restartOpt
.
Package neighbours (also on CRAN) offers helpers
for creating neighbourhood functions.
Examples
## Aim: find the columns of X that, when summed, give y
## random data set
nc < 25L ## number of columns in data set
nr < 5L ## number of rows in data set
howManyCols < 5L ## length of true solution
X < array(runif(nr*nc), dim = c(nr, nc))
xTRUE < sample(1L:nc, howManyCols)
Xt < X[ , xTRUE, drop = FALSE]
y < rowSums(Xt)
## a random solution x0 ...
makeRandomSol < function(nc) {
ii < sample.int(nc, sample.int(nc, 1L))
x0 < logical(nc); x0[ii] < TRUE
x0
}
x0 < makeRandomSol(nc)
## ... but probably not a good one
sum(y  rowSums(X[ , xTRUE, drop = FALSE])) ## should be 0
sum(y  rowSums(X[ , x0, drop = FALSE]))
## a neighbourhood function: switch n elements in solution
neighbour < function(xc, Data) {
xn < xc
p < sample.int(Data$nc, Data$n)
xn[p] < !xn[p]
if (sum(xn) < 1L)
xn < xc
xn
}
## a greedy neighbourhood function
neighbourG < function(xc, Data) {
of < function(x)
abs(sum(Data$y  rowSums(Data$X[ ,x, drop = FALSE])))
xbest < xc
Fxbest < of(xbest)
for (i in 1L:Data$nc) {
xn < xc; p < i
xn[p] < !xn[p]
if (sum(xn) >= 1L) {
Fxn < of(xn)
if (Fxn < Fxbest) {
xbest < xn
Fxbest < Fxn
}
}
}
xbest
}
## an objective function
OF < function(xn, Data)
abs(sum(Data$y  rowSums(Data$X[ ,xn, drop = FALSE])))
## (1) GREEDY SEARCH
## note: this could be done in a simpler fashion, but the
## redundancies/overhead here are small, and the example is to
## show how LSopt can be used for such a search
Data < list(X = X, y = y, nc = nc, nr = nr, n = 1L)
algo < list(nS = 500L, neighbour = neighbourG, x0 = x0,
printBar = FALSE, printDetail = FALSE)
solG < LSopt(OF, algo = algo, Data = Data)
## after how many iterations did we stop?
iterG < min(which(solG$Fmat[ ,2L] == solG$OFvalue))
solG$OFvalue ## the true solution has OFvalue 0
## (2) LOCAL SEARCH
algo$neighbour < neighbour
solLS < LSopt(OF, algo = algo, Data = Data)
iterLS < min(which(solLS$Fmat[ ,2L] == solLS$OFvalue))
solLS$OFvalue ## the true solution has OFvalue 0
## (3) *Threshold Accepting*
algo$nT < 10L
algo$nS < ceiling(algo$nS/algo$nT)
algo$q < 0.99
solTA < TAopt(OF, algo = algo, Data = Data)
iterTA < min(which(solTA$Fmat[ ,2L] == solTA$OFvalue))
solTA$OFvalue ## the true solution has OFvalue 0
## look at the solution
all < sort(unique(c(which(solTA$xbest),
which(solLS$xbest),
which(solG$xbest),
xTRUE)))
ta < ls < greedy < true < character(length(all))
true[ match(xTRUE, all)] < "o"
greedy[match(which(solG$xbest), all)] < "o"
ls[ match(which(solLS$xbest), all)] < "o"
ta[ match(which(solTA$xbest), all)] < "o"
data.frame(true = true, greedy = greedy, LS = ls , TA = ta,
row.names=all)
## plot results
par(ylog = TRUE, mar = c(5,5,1,6), las = 1)
plot(solTA$Fmat[seq_len(iterTA) ,2L],type = "l", log = "y",
ylim = c(1e4,
max(pretty(c(solG$Fmat,solLS$Fmat,solTA$Fmat)))),
xlab = "iterations", ylab = "OF value", col = grey(0.5))
lines(cummin(solTA$Fmat[seq_len(iterTA), 2L]), type = "l")
lines(solG$Fmat[ seq_len(iterG), 2L], type = "p", col = "blue")
lines(solLS$Fmat[seq_len(iterLS), 2L], type = "l", col = "goldenrod3")
legend(x = "bottomleft",
legend = c("TA best solution", "TA current solution",
"Greedy", "LS current/best solution"),
lty = c(1,1,0,1),
col = c("black",grey(0.5),"blue","goldenrod2"),
pch = c(NA,NA,21,NA))
axis(4, at = c(solG$OFvalue, solLS$OFvalue, solTA$OFvalue),
labels = NULL, las = 1)
lines(x = c(iterG, par()$usr[2L]), y = rep(solG$OFvalue,2),
col = "blue", lty = 3)
lines(x = c(iterTA, par()$usr[2L]), y = rep(solTA$OFvalue,2),
col = "black", lty = 3)
lines(x = c(iterLS, par()$usr[2L]), y = rep(solLS$OFvalue,2),
col = "goldenrod3", lty = 3)
Simple Moving Average
Description
The function computes a moving average of a vector.
Usage
MA(y, order, pad = NULL)
Arguments
y 
a numeric vector 
order 
An integer. The order of the moving average. The function is defined
such that order one returns 
pad 
Defaults to 
Value
Returns a vector of length length(y)
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
MA(1:10, 3)
MA(1:10, 3, pad = NA)
y < seq(1, 4, by = 0.3)
z < MA(y, 1)
all(y == z) ### (typically) FALSE
all.equal(y, z) ### should be TRUE
## 'Relative strength index'
rsi < function(y, t) {
y < diff(y)
ups < y + abs(y)
downs < y  abs(y)
RS < MA(ups, t) / MA(downs, t)
RS/(1 + RS)
}
x < cumprod(c(100, 1 + rnorm(100, sd = 0.01)))
par(mfrow = c(2,1))
plot(x, type = "l")
plot(rsi(x, 14), ylim = c(0,1), type = "l")
MaximumSharpeRatio/Tangency Portfolio
Description
Compute maximum Sharperatio portfolios, subject to lower and upper bounds on weights.
Usage
maxSharpe(m, var, min.return,
wmin = Inf, wmax = Inf, method = "qp",
groups = NULL, groups.wmin = NULL, groups.wmax = NULL)
Arguments
m 
vector of expected (excess) returns. 
var 
the covariance matrix: a numeric (real), symmetric matrix 
min.return 
minimumm required return. This is a technical parameter, used only for QP. 
wmin 
numeric: a lower bound on weights. May also be a vector that holds specific bounds for each asset. 
wmax 
numeric: an upper bound on weights. May also be a vector that holds specific bounds for each asset. 
method 
character. Currently, only 
groups 
a list of group definitions 
groups.wmin 
a numeric vector 
groups.wmax 
a numeric vector 
Details
The function uses solve.QP
from package
quadprog. Because of the algorithm that
solve.QP
uses, var
has to be positive
definit (i.e. must be of full rank).
Value
a numeric vector (the portfolio weights) with an attribute
variance
(the portfolio's variance)
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Schumann, E. (2012) Computing the global minimumvariance portfolio. http://enricoschumann.net/R/minvar.htm
See Also
minvar
,
mvPortfolio
,
mvFrontier
Examples
S < var(R < NMOF::randomReturns(3, 10, 0.03))
x < solve(S, colMeans(R))
x/sum(x)
x < coef(lm(rep(1, 10) ~ 1 + R))
unname(x/sum(x))
maxSharpe(m = colMeans(R), var = S)
maxSharpe(m = colMeans(R), var = S, wmin = 0, wmax = 1)
Option Pricing via MonteCarlo Simulation
Description
Functions to calculate the theoretical prices of options through simulation.
Usage
gbm(npaths, timesteps, r, v, tau, S0,
exp.result = TRUE, antithetic = FALSE)
gbb(npaths, timesteps, S0, ST, v, tau,
log = FALSE, exp.result = TRUE)
Arguments
npaths 
the number of paths 
timesteps 
timesteps per path 
r 
the mean per unit of time 
v 
the variance per unit of time 
tau 
time 
S0 
initial value 
ST 
final value of Brownian bridge 
log 
logical: construct bridge from log series? 
exp.result 
logical: compute 
antithetic 
logical: if 
Details
gbm
generates sample paths of geometric Brownian motion.
gbb
generates sample paths of a Brownian bridge by first creating
paths of Brownian motion W
from time 0
to time T
,
with W_0
equal to zero. Then, at each t
, it subtracts t/T
* W_T
and adds S0*(1t/T)+ST*(t/T)
.
Value
A matrix of sample paths; each column contains the price path of an
asset. Even with only a single timestep, the matrix will have two
rows (the first row is S0
).
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## price a European option
## ... parameters
npaths < 5000 ## increase number to get more precise results
timesteps < 1
S0 < 100
ST < 100
tau < 1
r < 0.01
v < 0.25^2
## ... create paths
paths < gbm(npaths, timesteps, r, v, tau, S0 = S0)
## ... a helper function
mc < function(paths, payoff, ...)
payoff(paths, ...)
## ... a payoff function (European call)
payoff < function(paths, X, r, tau)
exp(r * tau) * mean(pmax(paths[NROW(paths), ]  X, 0))
## ... compute and check
mc(paths, payoff, X = 100, r = r, tau = tau)
vanillaOptionEuropean(S0, X = 100, tau = tau, r = r, v = v)$value
## compute delta via forward difference
## (see Gilli/Maringer/Schumann, ch. 9)
h < 1e6 ## a small number
rnorm(1) ## make sure RNG is initialised
rnd.seed < .Random.seed ## store current seed
paths1 < gbm(npaths, timesteps, r, v, tau, S0 = S0)
.Random.seed < rnd.seed
paths2 < gbm(npaths, timesteps, r, v, tau, S0 = S0 + h)
delta.mc < (mc(paths2, payoff, X = 100, r = r, tau = tau)
mc(paths1, payoff, X = 100, r = r, tau = tau))/h
delta < vanillaOptionEuropean(S0, X = 100, tau = tau,
r = r, v = v)$delta
delta.mc  delta
## a fanplot
steps < 100
paths < results < gbm(1000, steps, r = 0, v = 0.2^2,
tau = 1, S0 = 100)
levels < seq(0.01, 0.49, length.out = 20)
greys < seq(0.9, 0.50, length.out = length(levels))
## start with an empty plot ...
plot(0:steps, rep(100, steps+1), ylim = range(paths),
xlab = "", ylab = "", lty = 0, type = "l")
## ... and add polygons
for (level in levels) {
l < apply(paths, 1, quantile, level)
u < apply(paths, 1, quantile, 1  level)
col < grey(greys[level == levels])
polygon(c(0:steps, steps:0), c(l, rev(u)),
col = col, border = NA)
## add border lines
## lines(0:steps, l, col = grey(0.4))
## lines(0:steps, u, col = grey(0.4))
}
Minimum ConditionalValueatRisk (CVaR) Portfolios
Description
Compute minimumCVaR portfolios, subject to lower and upper bounds on weights.
Usage
minCVaR(R, q = 0.1, wmin = 0, wmax = 1,
min.return = NULL, m = NULL,
method = "Rglpk",
groups = NULL, groups.wmin = NULL, groups.wmax = NULL,
Rglpk.control = list())
Arguments
R 
the scenario matrix: a numeric (real) matrix 
q 
the ValueatRisk level: a number between 0 and 0.5 
wmin 
numeric: a lower bound on weights. May also be a vector that holds specific bounds for each asset. 
wmax 
numeric: an upper bound on weights. May also be a vector that holds specific bounds for each asset. 
m 
vector of expected returns. Only used if 
min.return 
minimal required return. If 
method 
character. Currently, only 
groups 
a list of group definitions 
groups.wmin 
a numeric vector 
groups.wmax 
a numeric vector 
Rglpk.control 
a list: settings passed to 
Details
Compute the minimum CVaR portfolio for a given scenario set. The default method uses the formulation as a Linear Programme, as described in Rockafellar/Uryasev (2000).
The function uses Rglpk_solve_LP
from package
Rglpk.
Value
a numeric vector (the portfolio weights); attached is an
attribute whose name matches the method
name
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Rockafellar, R. T. and Uryasev, S. (2000). Optimization of Conditional ValueatRisk. Journal of Risk. 2 (3), 21–41.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Schumann, E. (2020) Minimising Conditional ValueatRisk (CVaR). http://enricoschumann.net/notes/minimisingconditionalvar.html
See Also
Examples
if (requireNamespace("Rglpk")) {
ns < 5000 ## number of scenarios
na < 20 ## nunber of assets
R < randomReturns(na, ns, sd = 0.01, rho = 0.5)
res < minCVaR(R, 0.25)
c(res) ## portfolio weights
}
Compute Minimum Mean–AbsoluteDeviation Portfolios
Description
Compute minimum mean–absolutedeviation portfolios.
Usage
minMAD(R, wmin = 0, wmax = 1,
min.return = NULL, m = NULL, demean = TRUE,
method = "lp",
groups = NULL, groups.wmin = NULL, groups.wmax = NULL,
Rglpk.control = list())
Arguments
R 
a matrix of return scenarios: each column represents one asset; each row represents one scenario 
wmin 
minimum weight 
wmax 
maximum weight 
min.return 
a minimum required return; ignored if 
m 
a vector of expected returns. If NULL, but 
demean 
logical. If 
method 
string. Supported are 
groups 
group definitions 
groups.wmin 
list of vectors 
groups.wmax 
list of vectors 
Rglpk.control 
a list 
Details
Compute the minimum mean–absolutedeviation portfolio for a given scenario set.
The function uses Rglpk_solve_LP
from package
Rglpk.
Value
a vector of portfolio weights
Author(s)
Enrico Schumann
References
Konno, H. and Yamazaki, H. (1991) MeanAbsolute Deviation Portfolio Optimization Model and Its Applications to Tokyo Stock Market. Management Science. 37 (5), 519–531.
See Also
Examples
na < 10
ns < 1000
R < randomReturns(na = na, ns = ns,
sd = 0.01, rho = 0.8, mean = 0.0005)
minMAD(R = R)
minvar(var(R))
MinimumVariance Portfolios
Description
Compute minimumvariance portfolios, subject to lower and upper bounds on weights.
Usage
minvar(var, wmin = 0, wmax = 1, method = "qp",
groups = NULL, groups.wmin = NULL, groups.wmax = NULL)
Arguments
var 
the covariance matrix: a numeric (real), symmetric matrix 
wmin 
numeric: a lower bound on weights. May also be a vector that holds specific bounds for each asset. 
wmax 
numeric: an upper bound on weights. May also be a vector that holds specific bounds for each asset. 
method 
character. Currently, only 
groups 
a list of group definitions 
groups.wmin 
a numeric vector 
groups.wmax 
a numeric vector 
Details
For method "qp"
, the function uses
solve.QP
from package
quadprog. Because of the algorithm that
solve.QP
uses, var
has to be positive
definite (i.e. must be of full rank).
Value
a numeric vector (the portfolio weights) with an attribute
variance
(the portfolio's variance)
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Schumann, E. (2012) Computing the global minimumvariance portfolio. http://enricoschumann.net/R/minvar.htm
See Also
Examples
## variancecovariance matrix from daily returns, 1 Jan 2014  31 Dec 2013, of
## cleaned data set at http://enricoschumann.net/data/gilli_accuracy.html
if (requireNamespace("quadprog")) {
var < structure(c(0.000988087100677907, 0.0000179669410403153, 0.000368923882626859,
0.000208303611101873, 0.000262742052359594, 0.0000179669410403153,
0.00171852167358765, 0.0000857467457561209, 0.0000215059246610556,
0.0000283532159921211, 0.000368923882626859, 0.0000857467457561209,
0.00075871953281751, 0.000194002299424151, 0.000188824454515841,
0.000208303611101873, 0.0000215059246610556, 0.000194002299424151,
0.000265780633005374, 0.000132611196599808, 0.000262742052359594,
0.0000283532159921211, 0.000188824454515841, 0.000132611196599808,
0.00025948420130626),
.Dim = c(5L, 5L),
.Dimnames = list(c("CBK.DE", "VOW.DE", "CON.DE", "LIN.DE", "MUV2.DE"),
c("CBK.DE", "VOW.DE", "CON.DE", "LIN.DE", "MUV2.DE")))
## CBK.DE VOW.DE CON.DE LIN.DE MUV2.DE
## CBK.DE 0.000988 0.0000180 0.0003689 0.0002083 0.0002627
## VOW.DE 0.000018 0.0017185 0.0000857 0.0000215 0.0000284
## CON.DE 0.000369 0.0000857 0.0007587 0.0001940 0.0001888
## LIN.DE 0.000208 0.0000215 0.0001940 0.0002658 0.0001326
## MUV2.DE 0.000263 0.0000284 0.0001888 0.0001326 0.0002595
##
minvar(var, wmin = 0, wmax = 0.5)
minvar(var,
wmin = c(0.1,0,0,0,0), ## enforce at least 10% weight in CBK.DE
wmax = 0.5)
minvar(var, wmin = Inf, wmax = Inf) ## no bounds
## [1] 0.0467 0.0900 0.0117 0.4534 0.4916
minvar(var, wmin = Inf, wmax = 0.45) ## no lower bounds
## [1] 0.0284 0.0977 0.0307 0.4500 0.4500
minvar(var, wmin = 0.1, wmax = Inf) ## no upper bounds
## [1] 0.100 0.100 0.100 0.363 0.337
## group constraints:
## group 1 consists of asset 1 only, and must have weight [0.25,0.30]
## group 2 consists of assets 4 and 5, and must have weight [0.10,0.20]
## => unconstrained
minvar(var, wmin = 0, wmax = 0.40)
## [1] 0.0097 0.1149 0.0754 0.4000 0.4000
## => with group constraints
minvar(var, wmin = 0, wmax = 0.40,
groups = list(1, 4:5),
groups.wmin = c(0.25, 0.1),
groups.wmax = c(0.30, 0.2))
## [1] 0.250 0.217 0.333 0.149 0.051
}
Computing Mean–Variance Efficient Portfolios
Description
Compute mean–variance efficient portfolios and efficient frontiers.
Usage
mvFrontier(m, var, wmin = 0, wmax = 1, n = 50, rf = NA,
groups = NULL, groups.wmin = NULL, groups.wmax = NULL)
mvPortfolio(m, var, min.return, wmin = 0, wmax = 1, lambda = NULL,
groups = NULL, groups.wmin = NULL, groups.wmax = NULL)
Arguments
m 
vector of expected returns 
var 
expected variance–covariance matrix 
wmin 
numeric: minimum weights 
wmax 
numeric: maximum weights 
n 
number of points on the efficient frontier 
min.return 
minimal required return 
rf 
riskfree rate 
lambda 
risk–reward tradeoff 
groups 
a list of group definitions 
groups.wmin 
a numeric vector 
groups.wmax 
a numeric vector 
Details
mvPortfolio
computes a single mean–variance
efficient portfolio, using package quadprog.
It does so by minimising portfolio variance, subject
to constraints on minimum return and budget (weights
need to sum to one), and min/max constraints on the
weights.
If \lambda
is specified, the function ignores the min.return
constraint and instead solves the model
\min_w\ \ \lambda \mbox{\code{m}}'w + (1\lambda)
w'\mbox{\code{var}\,}w
in which w
are the weights. If
\lambda
is a vector of length 2, then the model becomes
\min_w\ \ \lambda_1 \mbox{\code{m}\,}'w + \lambda_2
w'\mbox{\code{var}\,}w
which may be more convenient
(e.g. for setting \lambda_1
to 1).
mvFrontier
computes returns, volatilities and
compositions for portfolios along an efficient frontier.
If rf
is not NA
, cash is included as an asset.
Value
For mvPortfolio
, a numeric vector of weights.
For mvFrontier
, a list of three components:
return 
returns of portfolios 
volatility 
volatilities of portfolios 
weights 
A matrix of portfolio weights.
Each column holds the weights for one portfolio on the
frontier. If 
The ith portfolio on the frontier corresponds
to the ith elements of return
and
volatility
, and the ith column of
portfolio
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
minvar
for computing the minimumvariance portfolio
Examples
na < 4
vols < c(0.10, 0.15, 0.20,0.22)
m < c(0.06, 0.12, 0.09, 0.07)
const_cor < function(rho, na) {
C < array(rho, dim = c(na, na))
diag(C) < 1
C
}
var < diag(vols) %*% const_cor(0.5, na) %*% diag(vols)
wmax < 1 # maximum holding size
wmin < 0.0 # minimum holding size
rf < 0.02
if (requireNamespace("quadprog")) {
p1 < mvFrontier(m, var, wmin = wmin, wmax = wmax, n = 50)
p2 < mvFrontier(m, var, wmin = wmin, wmax = wmax, n = 50, rf = rf)
plot(p1$volatility, p1$return, pch = 19, cex = 0.5, type = "o",
xlab = "Expected volatility",
ylab = "Expected return")
lines(p2$volatility, p2$return, col = grey(0.5))
abline(v = 0, h = rf)
} else
message("Package 'quadprog' is required")
Internal NMOF functions
Description
Several internal functions, not exported from the NMOF name space.
Usage
makeInteger(x, label, min = 1L)
anyNA(x)
checkList(passedList, defaultList, label = "'algo'")
mRU(m,n)
mRN(m,n)
mcList(mc.control)
repair1c(x, up, lo)
Chapters1
Chapters2
due(D, tauD, tau, q)
Details
makeInteger
coerces the scalar x
to
integer. If the result is NA
or smaller than
min
, an error is issued.
anyNA
gives TRUE
if any(is.na(x))
, else
FALSE
. If x
is a function or NULL
, it
also gives FALSE
.
checkList
will issue an error if passedList
contains any NA
elements. It will give a warning if
any elements in passedList
are unnamed, or if an
element in names(passedList)
is not found in
names(defaultList)
.
mRU
and mRN
create matrices (of size m
times n
) of uniform/Gaussian variates.
mcList
takes a list of named elements and
‘merges’ them with the default settings of
mclapply
from package parallel.
repair1c
is described in the vignette on
‘Repairing Solutions’.
Chapters[12]
is a vector of length 15, giving the
chapter names as printed in the book. (Can be accessed with
showChapterNames
.)
due
(dividends until expiry) returns a list with
named components tauD
and D
: all dividends for
which timetopayment is greater than zero, but less than or
equal to timetoexpiry.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Zero Rates for Nelson–Siegel–Svensson Model
Description
Compute zero yields for Nelson–Siegel (NS)/Nelson–Siegel–Svensson (NSS) model.
Usage
NS(param, tm)
NSS(param, tm)
Arguments
param 
a vector. For NS: 
tm 
a vector of maturities 
Details
See Chapter 14 in Gilli/Maringer/Schumann (2011).
Maturities (tm
) need to be given in time (not dates).
Value
The function returns a vector of length length(tm)
.
Author(s)
Enrico Schumann
References
Gilli, M. and Grosse, S. and Schumann, E. (2010) Calibrating the NelsonSiegelSvensson model, COMISEF Working Paper Series No. 031. http://enricoschumann.net/COMISEF/wps031.pdf
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Gilli, M. and Schumann, E. (2010) A Note on ‘Good’ Starting Values in Numerical Optimisation, COMISEF Working Paper Series No. 044. http://enricoschumann.net/COMISEF/wps044.pdf
Nelson, C.R. and Siegel, A.F. (1987) Parsimonious Modeling of Yield Curves. Journal of Business, 60(4), pp. 473–489.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Svensson, L.E. (1994) Estimating and Interpreting Forward Interest Rates: Sweden 1992–1994. IMF Working Paper 94/114.
See Also
Examples
tm < c(c(1, 3, 6, 9) / 12, 1:10) ## in years
param < c(6, 3, 8, 1)
yM < NS(param, tm)
plot(tm, yM, xlab = "maturity in years",
ylab = "yield in percent")
param < c(6, 3, 5, 5, 1, 3)
yM < NSS(param, tm)
plot(tm, yM, xlab = "maturity in years",
ylab = "yield in percent")
## get Bliss/Diebold/Li data (used in some of the papers in References)
u < url("https://www.sas.upenn.edu/~fdiebold/papers/paper49/FBFITTED.txt")
try(open(u))
BliDiLi < try(scan(u, skip = 14))
if (!inherits(BliDiLi, "tryerror")) {
close(u)
mat < NULL
for (i in 1:372)
mat < rbind(mat,BliDiLi[(19*(i1)+1):(19*(i1)+19)])
mats < c(1,3,6,9,12,15,18,21,24,30,36,48,60,72,84,96,108,120)/12
## the obligatory perspective plot
persp(x = mat[,1], y = mats, mat[ ,1L],
phi = 30, theta = 30, ticktype = "detailed",
xlab = "time",
ylab = "time to maturity in years",
zlab = "zero rates in %")
}
Factor Loadings for Nelson–Siegel and Nelson–Siegel–Svensson
Description
Computes the factor loadings for Nelson–Siegel (NS) and Nelson–Siegel–Svensson (NSS) model for given lambda
values.
Usage
NSf(lambda, tm)
NSSf(lambda1, lambda2, tm)
Arguments
lambda 
the 
lambda1 
the 
lambda2 
the 
tm 
a numeric vector with timestopayment/maturity 
Details
The function computes the factor loadings for given \lambda
parameters. Checking the correlation between these
factor loadings can help to set reasonable \lambda
values for the NS/NSS models.
Value
For NS, a matrix with length(tm)
rows and three columns.
For NSS, a matrix with length(tm)
rows and four columns.
Author(s)
Enrico Schumann
References
Gilli, M. and Grosse, S. and Schumann, E. (2010) Calibrating the NelsonSiegelSvensson model, COMISEF Working Paper Series No. 031. http://enricoschumann.net/COMISEF/wps031.pdf
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Gilli, M. and Schumann, E. (2010) A Note on ‘Good’ Starting Values in Numerical Optimisation, COMISEF Working Paper Series No. 044. http://enricoschumann.net/COMISEF/wps044.pdf
Nelson, C.R. and Siegel, A.F. (1987) Parsimonious Modeling of Yield Curves. Journal of Business, 60(4), pp. 473–489.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Svensson, L.E. (1994) Estimating and Interpreting Forward Interest Rates: Sweden 1992–1994. IMF Working Paper 94/114.
See Also
Examples
## NelsonSiegel
cor(NSf(lambda = 6, tm = 1:10)[1L, 1L])
## NelsonSiegelSvensson
cor(NSSf(lambda1 = 1, lambda2 = 5, tm = 1:10)[1L, 1L])
cor(NSSf(lambda1 = 4, lambda2 = 9, tm = 1:10)[1L, 1L])
Option Data
Description
Closing prices of DAX index options as of 20120210.
Usage
optionData
Format
optionData
is a list with six components:
pricesCall
a matrix of size 124 times 10. The rows are the strikes; each column belongs to one expiry date.
pricesPut
a matrix of size 124 times 10
index
The DAX index (spot).
future
The available future settlement prices.
Euribor
Euribor rates.
NSSpar
Paramaters for German government bond yields, as estimated by the Bundesbank.
Details
Settlement prices for EUREX options are computed at 17:30, Frankfurt Time, even though trading continues until 22:00.
Source
The data was obtained from several websites: close prices of EUREX products were collected from https://www.eurex.com/exen/ ; Euribor rates and the parameters of the NelsonSiegelSvensson can be found at https://www.bundesbank.de/en/ .
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
str(optionData)
NSS(optionData$NSSpar, 1:10)
Partial Moments
Description
Compute partial moments.
Usage
pm(x, xp = 2, threshold = 0, lower = TRUE,
normalise = FALSE, na.rm = FALSE)
Arguments
x 
a numeric vector or a matrix 
xp 
exponent 
threshold 
a numeric vector of length one 
lower 
logical 
normalise 
logical 
na.rm 
logical 
Details
For a vector x
of length n
, partial
moments are computed as follows:
\mathrm{upper\ partial\ moment} = \frac{1}{n} \sum_{x >
t}\left(x  t \right)^e
\mathrm{lower\ partial\ moment} = \frac{1}{n} \sum_{x <
t}\left(t  x \right)^e
The threshold
is denoted t
, the exponent
xp
is labelled e
.
If normalise
is TRUE
, the result is raised to
1/xp
. If x
is a matrix, the function will compute the
partial moments columnwise.
See Gilli, Maringer and Schumann (2019), chapter 14.
Value
numeric
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
pm(x < rnorm(100), 2)
var(x)/2
pm(x, 2, normalise = TRUE)
sqrt(var(x)/2)
Particle Swarm Optimisation
Description
The function implements Particle Swarm Optimisation.
Usage
PSopt(OF, algo = list(), ...)
Arguments
OF 
the objective function to be minimised. See Details. 
algo 
a list with the settings for algorithm. See Details and Examples. 
... 
pieces of data required to evaluate the objective function. See Details. 
Details
The function implements Particle Swarm Optimisation (PS); see the references for details on the implementation. PS is a populationbased optimisation heuristic. It develops several solutions (a ‘population’) over a number of iterations. PS is directly applicable to continuous problems since the population is stored in realvalued vectors. In each iteration, a solution is updated by adding another vector called velocity. Think of a solution as a position in the search space, and of velocity as the direction into which this solution moves. Velocity changes over the course of the optimization: it is biased towards the best solution found by the particular solution and the best overall solution. The algorithm stops after a fixed number of iterations.
To allow for constraints, the evaluation works as follows: after a new
solution is created, it is (i) repaired, (ii) evaluated through the
objective function, (iii) penalised. Step (ii) is done by a call to
OF
; steps (i) and (iii) by calls to algo$repair
and
algo$pen
. Step (i) and (iii) are optional, so the respective
functions default to NULL
. A penalty can also be directly
written in the OF
, since it amounts to a positive number added
to the ‘clean’ objective function value. It can be
advantageous to write a separate penalty function if either only the
objective function or only the penalty function can be vectorised.
(Constraints can also be added without these mechanisms. Solutions
that violate constraints can, for instance, be mapped to feasible
solutions, but without actually changing them. See Maringer and
Oyewumi, 2007, for an example with Differential Evolution.)
Conceptually, PS consists of two loops: one loop across the
iterations and, in any given generation, one loop across the
solutions. This is the default, controlled by the variables
algo$loopOF
, algo$loopRepair
, algo$loopPen
and
loopChangeV
which all default to TRUE
. But it does not
matter in what order the solutions are evaluated (or repaired or
penalised), so the second loop can be vectorised. Examples are given
in the vignettes and in the book. The respective algo$loopFun
must then be set to FALSE
.
The objective function, the repair function and and the penalty
function will be called as fun(solution, ...)
.
The list algo
contains the following items:
nP
population size. Defaults to 100. Using default settings may not be a good idea.
nG
number of iterations. Defaults to 500. Using default settings may not be a good idea.
c1
the weight towards the individual's best solution. Typically between 0 and 2; defaults to 1. Using default settings may not be a good idea. In some cases, even negative values work well: the solution is then driven off its past best position. For ‘simple’ problems, setting
c1
to zero may work well: the population moves then towards the best overall solution.c2
the weight towards the populations's best solution. Typically between 0 and 2; defaults to 1. Using default settings may not be a good idea. In some cases, even negative values work well: the solution is then driven off the population's past best position.
iner
the inertia weight (a scalar), which reduces velocity. Typically between 0 and 1. Default is 0.9.
initV
the standard deviation of the initial velocities. Defaults to 1.
maxV
the maximum (absolute) velocity. Setting limits to velocity is sometimes called velocity clamping. Velocity is the change in a given solution in a given iteration. A maximum velocity can be set so to prevent unreasonable velocities (‘overshooting’): for instance, if a decision variable may lie between 0 and 1, then an absolute velocity much greater than 1 makes rarely sense.
min
,max
vectors of minimum and maximum parameter values. The vectors
min
andmax
are used to determine the dimension of the problem and to randomly initialise the population. Per default, they are no constraints: a solution may well be outside these limits. Only ifalgo$minmaxConstr
isTRUE
will the algorithm repair solutions outside themin
andmax
range.minmaxConstr
if
TRUE
,algo$min
andalgo$max
are considered constraints. Default isFALSE
.pen
a penalty function. Default is
NULL
(no penalty).repair
a repair function. Default is
NULL
(no repairing).changeV
a function to change velocity. Default is
NULL
(no change). This function is called before the velocity is added to the current solutions; it can be used to impose restrictions like changing only a number of decision variables.initP
optional: the initial population. A matrix of size
length(algo$min)
timesalgo$nP
, or a function that creates such a matrix. If a function, it should take no arguments.loopOF
logical. Should the
OF
be evaluated through a loop? Defaults toTRUE
.loopPen
logical. Should the penalty function (if specified) be evaluated through a loop? Defaults to
TRUE
.loopRepair
logical. Should the repair function (if specified) be evaluated through a loop? Defaults to
TRUE
.loopChangeV
logical. Should the
changeV
function (if specified) be evaluated through a loop? Defaults toTRUE
.printDetail
If
TRUE
(the default), information is printed. If an integeri
greater then one, information is printed at veryi
th iteration.printBar
If
TRUE
(the default), atxtProgressBar
(from package utils) is printed).storeF
If
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
default is
FALSE
. IfTRUE
, the solutions (ie, decision variables) in every generation are stored as listsP
andPbest
, both stored in the listxlist
which the function returns. To check, for instance, the solutions at the end of thei
th iteration, retrievexlist[[c(1L, i)]]
; the best solutions at the end of this iteration are inxlist[[c(2L, i)]]
.P[[i]]
andPbest[[i]]
will be matrices of sizelength(algo$min)
timesalgo$nP
.classify
Logical; default is
FALSE
. IfTRUE
, the result will have a class attributeTAopt
attached. This feature is experimental: the supported methods may change without warning.drop

Default is
TRUE
. IfFALSE
, the dimension is not dropped from a single solution when it is passed to a function. (That is, the function will receive a singlecolumn matrix.)
Value
Returns a list:
xbest 
the solution 
OFvalue 
objective function value of best solution 
popF 
a vector: the objective function values in the final population 
Fmat 
if 
xlist 
if 
initial.state 
the value of 
Author(s)
Enrico Schumann
References
Eberhart, R.C. and Kennedy, J. (1995) A New Optimizer using Particle Swarm theory. Proceedings of the Sixth International Symposium on Micromachine and Human Science, pp. 39–43.
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## Least Median of Squares (LMS) estimation
genData < function(nP, nO, ol, dy) {
## create dataset as in SalibianBarrera & Yohai 2006
## nP = regressors, nO = number of obs
## ol = number of outliers, dy = outlier size
mRN < function(m, n) array(rnorm(m * n), dim = c(m, n))
y < mRN(nO, 1)
X < cbind(as.matrix(numeric(nO) + 1), mRN(nO, nP  1L))
zz < sample(nO)
z < cbind(1, 100, array(0, dim = c(1L, nP  2L)))
for (i in seq_len(ol)) {
X[zz[i], ] < z
y[zz[i]] < dy
}
list(X = X, y = y)
}
OF < function(param, data) {
X < data$X
y < data$y
aux < as.vector(y)  X %*% param
## as.vector(y) for recycling (param is a matrix)
aux < aux * aux
aux < apply(aux, 2, sort, partial = data$h)
aux[h, ]
}
nP < 2L; nO < 100L; ol < 10L; dy < 150
aux < genData(nP,nO,ol,dy); X < aux$X; y < aux$y
h < (nO + nP + 1L) %/% 2
data < list(y = y, X = X, h = h)
algo < list(min = rep(10, nP), max = rep( 10, nP),
c1 = 1.0, c2 = 2.0,
iner = 0.7, initV = 1, maxV = 3,
nP = 100L, nG = 300L, loopOF = FALSE)
system.time(sol < PSopt(OF = OF, algo = algo, data = data))
if (require("MASS", quietly = TRUE)) {
## for nsamp = "best", in this case, complete enumeration
## will be tried. See ?lqs
system.time(test1 < lqs(data$y ~ data$X[, 1L],
adjust = TRUE,
nsamp = "best",
method = "lqs",
quantile = data$h))
}
## check
x1 < sort((y  X %*% as.matrix(sol$xbest))^2)[h]
cat("Particle Swarm\n",x1,"\n\n")
if (require("MASS", quietly = TRUE)) {
x2 < sort((y  X %*% as.matrix(coef(test1)))^2)[h]
cat("lqs\n", x2, "\n\n")
}
PutCall Parity
Description
Put–call parity
Usage
putCallParity(what, call, put, S, X, tau, r, q = 0, tauD = 0, D = 0)
Arguments
what 
character: what to compute. Currently only 
call 
call price 
put 
put price 
S 
underlier 
X 
strike 
tau 
time to expiry 
r 
interest rate 
q 
dividend rate 
tauD 
numeric vector: time to dividend 
D 
numeric vector: dividends 
Details
Put–call parity only works for European options. The function is
vectorised (like vanillaOptionEuropean
), except for
dividends.
Value
Numeric vector.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
S < 100; X < 100; tau < 1; r < 0.02; q < 0.0;
vol < 0.3; D < 20; tauD < 0.5
call < vanillaOptionEuropean(S, X, tau, r, q, vol^2,
tauD = tauD, D = D, type = "call")$value
put < vanillaOptionEuropean(S, X, tau, r, q, vol^2,
tauD = tauD, D = D, type = "put")$value
## recover the call from the put (et vice versa)
all.equal(call, putCallParity("call", put = put, S=S, X=X, tau=tau,
r=r, q=q, tauD=tauD, D=D))
all.equal(put, putCallParity("put", call = call, S=S, X=X, tau=tau,
r=r, q=q, tauD=tauD, D=D))
## BlackScholesMerton with with 'callCF'
S < 100; X < 90; tau < 1; r < 0.02; q < 0.08
v < 0.2^2 ## variance, not volatility
(ccf < callCF(cf = cfBSM, S = S, X = X, tau = tau, r = r, q = q,
v = v, implVol = TRUE))
all.equal(ccf$value,
vanillaOptionEuropean(S, X, tau, r, q, v, type = "call")$value)
all.equal(
putCallParity("put", call=ccf$value, S=S, X=X, tau=tau, r=r, q=q),
vanillaOptionEuropean(S, X, tau, r, q, v, type = "put")$value)
Prepare LaTeX Table with Quartile Plots
Description
The function returns the skeleton of a LaTeX tabular that contains the
median, minimum and maximum of the columns of a matrix X
. For
each column, a quartile plot is added.
Usage
qTable(X, xmin = NULL, xmax = NULL, labels = NULL, at = NULL,
unitlength = "5cm", linethickness = NULL,
cnames = colnames(X), circlesize = 0.01,
xoffset = 0, yoffset = 0, dec = 2, filename = NULL,
funs = list(median = median, min = min, max = max),
tabular.format, skip = TRUE)
Arguments
X 
a numeric matrix (or an object that can be coerced to a numeric
matrix with 
xmin 
optional: the minimum for the xaxis. See Details. 
xmax 
optional: the maximum for the xaxis. See Details. 
labels 
optional: labels for the xaxis. 
at 
optional: where to put labels. 
unitlength 
the unitlength for LaTeX's 
linethickness 
the linethickness for LaTeX's 
cnames 
the column names of 
circlesize 
the size of the circle in LaTeX's 
xoffset 
defaults to 0. See Details. 
yoffset 
defaults to 0. See Details. 
dec 
the number of decimals 
filename 
if provided, output is 
funs 
A

tabular.format 
optional: character string like 
skip 
Adds a newline at the end of the tabular. Default is

Details
The function creates a onecolumn character matrix that can be put into
a LaTeX file (the matrix holds a tabular). It relies on LaTeX's
picture
environment and should work for LaTeX and pdfLaTeX. Note
that the tabular needs generally be refined, depending on the settings
and the data.
The tabular has one row for every column of X
(and header and
footer rows). A given row contains (per default) the median, the minimum
and the maximum of the column; it also includes a picture
environment the shows a quartile plot of the distribution of the
elements in that column. Other functions can be specified via argument
funs
.
A number of parameters can be passed to LaTeX's picture
environment: unitlength
, xoffset
, yoffset
,
linethickness
. Sizes and lengths are functions of
unitlength
(linethickness
is an exception; and while
circlesize
is a multiple of unitlength, it will not translate
into an actual diameter of more than 14mm).
The whole tabular environment is put into curly brackets so that the settings do not change settings elsewhere in the LaTeX document.
If xmin
, xmax
, labels
and at
are not
specified, they are computed through a call to pretty
from
the base package. If limits are specified, then both xmin
and xmax
must be set; if labels are used, then both labels
and at
must be specified.
To use the function in a vignette, use cat(tTable(X))
(and
results=tex
in the code chunk options). The vignette
qTableEx
shows some examples.
Value
A matrix of mode character. If filename
is specified then
qTable
will have the side effect of writing a textfile with a
LaTeX tabular.
Note
qTable
returns a raw draft of a table for LaTeX. Please, spend
some time on making it pretty.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Tufte, E. (2001) The Visual Display of Quantitative Information. 2nd edition, Graphics Press.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
x < rnorm(100, mean = 0, sd = 2)
y < rnorm(100, mean = 1, sd = 2)
z < rnorm(100, mean = 1, sd = 0.5)
X < cbind(x, y, z)
res < qTable(X)
print(res)
cat(res)
## Not run:
## show vignette with examples
qt < vignette("qTableEx", package = "NMOF")
print(qt)
edit(qt)
## create a simple LaTeX file 'test.tex':
## 
## \documentclass{article}
## \begin{document}
## \input{res.tex}
## \end{document}
## 
res < qTable(X, filename = "res.tex", yoffset = 0.025, unitlength = "5cm",
circlesize = 0.0125, xmin = 10, xmax = 10, dec = 2)
## End(Not run)
Create a Random Returns
Description
Create a matrix of random returns.
Usage
randomReturns(na, ns, sd, mean = 0, rho = 0, exact = FALSE)
Arguments
na 
number of assets 
ns 
number of return scenarios 
sd 
the standard deviation: either a single number or a vector
of length 
mean 
the mean return: either a single number or a vector
of length 
rho 
correlation: either a scalar (i.e. a constant pairwise correlation) or a correlation matrix 
exact 
logical: if 
Details
The function corresponds to the function random_returns
,
described in the second edition of NMOF (the book).
Value
a numeric
matrix
of size na
times
ns
Note
The function corresponds to the function random_returns
,
described in the second edition of NMOF (the book).
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
if (requireNamespace("quadprog")) {
## a small experiment: when computing minimumvariance portfolios
## for correlated assets, how many large positions are in the portfolio?
na < 100 ## number of assets
inc < 5 ## minimum of assets to include
n < numeric(10)
for (i in seq_along(n)) {
R < randomReturns(na = na,
ns = 500,
sd = seq(.2/.16, .5/.16, length.out = 100),
rho = 0.5)
n[i] < sum(minvar(cov(R), wmax = 1/inc)> 0.01)
}
summary(n)
}
Repair an Indefinite Correlation Matrix
Description
The function ‘repairs’ an indefinite correlation matrix by replacing its negative eigenvalues by zero.
Usage
repairMatrix(C, eps = 0)
Arguments
C 
a correlation matrix 
eps 
a small number 
Details
The function ‘repairs’ a correlation matrix: it
replaces negative eigenvalues with eps
and rescales
the matrix such that all elements on the main diagonal
become unity again.
Value
Returns a numeric matrix.
Note
This function may help to cure a numerical problem, but it will rarely help to cure an empirical problem. (Garbage in, garbage out.)
See also the function nearPD
in the Matrix package.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Rebonato, R. and Jaeckel, P. (1999) The most general methodology to create a valid correlation matrix for risk management and option pricing purposes.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
## example: build a portfolio of three assets
C < c(1,.9,.9,.9,1,.2,.9,.2,1)
dim(C) < c(3L, 3L)
eigen(C, only.values = TRUE)
vols < c(.3, .3, .3) ## volatilities
S < C * outer(vols,vols) ## covariance matrix
w < c(1, 1, 1) ## a portfolio
w %*% S %*% w ## variance of portfolio is negative!
sqrt(as.complex(w %*% S %*% w))
S < repairMatrix(C) * outer(vols,vols)
w %*% S %*% w ## more reasonable
sqrt(w %*% S %*% w)
Resample with Specified Rank Correlation
Description
Resample with replacement from a number of vectors; the sample will have a specified rank correlation.
Usage
resampleC(..., size, cormat)
Arguments
... 
numeric vectors; they need not have the same length. 
size 
an integer: the number of samples to draw 
cormat 
the rank correlation matrix 
Details
See Gilli, Maringer and Schumann (2011), Section 7.1.2. The function
samples with replacement from the vectors passed through
...
. The resulting samples will have an (approximate) rank
correlation as specified in cormat
.
The function uses the eigenvalue decomposition to generate the
correlation; it will not break down in case of a semidefinite
matrix. If an eigenvalue of cormat
is smaller than zero, a
warning is issued (but the function proceeds).
Value
a numeric matrix with size
rows. The columns contain the
samples; hence, there will be as many columns as vectors passed
through ...
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## a sample
v1 < rnorm(20)
v2 < runif(50)
v3 < rbinom(100, size = 50, prob = 0.4)
## a correlation matrix
cormat < array(0.5, dim = c(3, 3))
diag(cormat) < 1
cor(resampleC(a = v1, b = v2, v3, size = 100, cormat = cormat),
method = "spearman")
Restart an Optimisation Algorithm
Description
The function provides a simple wrapper for the optimisation algorithms in the package.
Usage
restartOpt(fun, n, OF, algo, ...,
method = c("loop", "multicore", "snow"),
mc.control = list(), cl = NULL,
best.only = FALSE)
Arguments
fun 
the optimisation function: 
n 
the number of restarts 
OF 
the objective function 
algo 
the list 
... 
additional data that is passed to the particular optimisation function 
method 
can be 
mc.control 
a list containing settings that will be passed to 
cl 
default is 
best.only 
if 
Details
The function returns a list of lists. If a specific starting solution
is passed, all runs will start from this solution. If this is not
desired, initial solutions can be created randomly. This is done per
default in DEopt
, GAopt
and
PSopt
, but LSopt
and TAopt
require to specify a starting solution.
In case of LSopt
and TAopt
, the passed
initial solution algo$x0
is checked with is.function
: if
TRUE
, the function is evaluated in each single run. For
DEopt
, GAopt
and PSopt
, the
initial solution (which also can be a function) is specified with
algo$initP
.
The argument method
determines how fun
is
evaluated. Default is loop
. If method
is "multicore",
function mclapply
from package parallel is used. Further
settings for mclapply
can be passed through the list
mc.control
. If multicore
is chosen but the functionality
is not available, then method
will be set to loop
and a
warning is issued. If method == "snow"
, function
clusterApply
from package parallel is used. In this case,
the argument cl
must either be a cluster object (see the
documentation of clusterApply
) or an integer. If an integer, a
cluster will be set up via makeCluster(c(rep("localhost", cl)),
type = "SOCK")
, and stopCluster
is called when the function is
exited. If snow
is chosen but parallel is not available
or cl
is not specified, then method
will be set to
loop
and a warning is issued. In case that cl
is an
cluster object, stopCluster
will not be called automatically.
Value
If best.only
is FALSE
(the default), the function
returns a list of n
lists. Each of the n
lists stores
the output of one of the runs.
If best.only
is TRUE
, only the best restart is
reported. The returned list has the structure specific to the used
method.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
DEopt
, GAopt
, LSopt
,
PSopt
, TAopt
Examples
## see example(DEopt)
algo < list(nP = 50L,
F = 0.5,
CR = 0.9,
min = c(10, 10),
max = c( 10, 10),
printDetail = FALSE,
printBar = FALSE)
## choose a larger 'n' when you can afford it
algo$nG < 100L
res100 < restartOpt(DEopt, n = 5L, OF = tfTrefethen, algo = algo)
res100F < sapply(res100, `[[`, "OFvalue")
algo$nG < 200L
res200 < restartOpt(DEopt, n = 5L, OF = tfTrefethen, algo = algo)
res200F < sapply(res200, `[[`, "OFvalue")
xx < pretty(c(res100F, res200F, 3.31))
plot(ecdf(res100F), main = "optimum is 3.306",
xlim = c(xx[1L], tail(xx, 1L)))
abline(v = 3.3069, col = "red") ## optimum
lines(ecdf(res200F), col = "blue")
legend(x = "right", box.lty = 0, , lty = 1,
legend = c("optimum", "100 generations", "200 generations"),
pch = c(NA, 19, 19), col = c("red", "black", "blue"))
## a 'bestofN' strategy: given a sample x of objective
## function values, compute the probability that, after N draws,
## we have at least one realisation not worse than X
x < c(0.1,.3,.5,.5,.6)
bestofN < function(x, N) {
nx < length(x)
function(X)
1  (sum(x > X)/nx)^N
}
bestof2 < bestofN(x, 2)
bestof5 < bestofN(x, 5)
bestof2(0.15)
bestof5(0.15)
## Not run:
## with R >= 2.13.0 and the compiler package
algo$nG < 100L
system.time(res100 < restartOpt(DEopt, n = 10L, OF = tfTrefethen, algo = algo))
require("compiler")
enableJIT(3)
system.time(res100 < restartOpt(DEopt, n = 10L, OF = tfTrefethen, algo = algo))
## End(Not run)
Download Jay Ritter's IPO Data
Description
Download IPO data provided by Jay Ritter and transform them into a data frame.
Usage
Ritter(dest.dir,
url = "https://site.warrington.ufl.edu/ritter/files/IPOage.xlsx")
Arguments
dest.dir 
character: a path to a directory 
url 
the data URL 
Details
The function downloads IPO data provided by Jay R. Ritter https://site.warrington.ufl.edu/ritter. Since the data are provided in Excel format, package openxlsx is required.
The downloaded Excel gets a date prefix (today in
format YYYYMMDD
) and is stored in directory
dest.dir
. Before any download is attempted,
the function checks whether a file with today's
prefix exist in dest.dir
; if yes, this file is
used.
Value
a data.frame
:
CUSIP 
CUSIP 
Offer date 
a 
Company name 
character: Company name 
Ticker 
character: Ticker 
Founding 
Founding year 
PERM 
PERM 
VC dummy 
VC Dummy 
Rollup 
Rollup 
Dual 
Dual 
Postissue shares 
Postissue shares 
Internet 
Internet 
Author(s)
Enrico Schumann
References
https://site.warrington.ufl.edu/ritter/ipodata/
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## Not run:
archive.dir < "~/Downloads/Ritter"
if (!dir.exists(archive.dir))
dir.create(archive.dir)
Ritter(archive.dir)
## End(Not run)
SimulatedAnnealing Information
Description
The function can be called from the objective and neighbourhood
function during a run of SAopt
; it provides information
such as the current iteration, the current solution, etc.
Usage
SA.info(n = 0L)
Arguments
n 
generational offset; see Details. 
Details
This function is still experimental.
The function can be called in the neighbourhood function or the
objective function during a run of SAopt
. It evaluates
to a list with information about the state of the optimisation run,
such as the current iteration or the currently best solution.
SA.info
relies on parent.frame
to retrieve its
information. If the function is called within another function within
the neighbourhood or objective function, the argument n
needs
to be increased.
Value
A list
calibration 
logical: whether the algorithm is calibrating the acceptance probability 
iteration 
current iteration 
step 
current step for the given temperature level 
temperature 
current temperature (the number, not the value) 
xbest 
the best solution found so far 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
### MINIMAL EXAMPLE for SAopt
## the objective function evaluates to a constant
fun < function(x)
0
## the neighbourhood function does not even change
## the solution; it only reports information
nb < function(x) {
info < SA.info()
cat("current step ", info$step,
" current iteration ", info$iteration, "\n")
x
}
## run SA
algo < list(nS = 5, nT = 2, nD = 10,
initT = 1,
x0 = rep(0, 5),
neighbour = nb,
printBar = FALSE)
ignore < SAopt(fun, algo)
Optimisation with Simulated Annealing
Description
The function implements a SimulatedAnnealing algorithm.
Usage
SAopt(OF, algo = list(), ...)
Arguments
OF 
The objective function, to be minimised. Its first argument
needs to be a solution 
algo 
A list of settings for the algorithm. See Details. 
... 
other variables passed to 
Details
Simulated Annealing (SA) changes an initial solution
iteratively; the algorithm stops after a fixed number of
iterations. Conceptually, SA consists of a loop than runs
for a number of iterations. In each iteration, a current solution
xc
is changed through a function algo$neighbour
. If this
new (or neighbour) solution xn
is not worse than xc
, ie,
if OF(xn,...) <= OF(xc,...)
, then xn
replaces
xc
. If xn
is worse, it still replaces xc
, but
only with a certain probability. This probability is a function of the
degree of the deterioration (the greater, the less likely the new
solution is accepted) and the current iteration (the longer the
algorithm has already run, the less likely the new
solution is accepted).
The list algo
contains the following items.
nS
The number of steps per temperature. The default is 1000; but this setting depends very much on the problem.
nT
The number of temperatures. Default is 10.
nI

Total number of iterations, with default
NULL
. If specified, it will overridenS
withceiling(nI/nT)
. Using this option makes it easier to compare and switch between functionsLSopt
,TAopt
andSAopt
. nD
The number of random steps to calibrate the temperature. Defaults to 2000.
initT
Initial temperature. Defaults to
NULL
, in which case it is automatically chosen so thatinitProb
is achieved.finalT
Final temperature. Defaults to 0.
alpha
The cooling constant. The current temperature is multiplied by this value. Default is 0.9.
mStep
Step multiplier. The default is 1, which implies constant number of steps per temperature. If greater than 1, the step number
nS
is increased tom*nS
(and rounded).x0
The initial solution. If this is a function, it will be called once without arguments to compute an initial solution, ie,
x0 < algo$x0()
. This can be useful when the routine is called in a loop of restarts, and each restart is to have its own starting value.neighbour
The neighbourhood function, called as
neighbour(x, ...)
. Its first argument must be a solutionx
; it must return a changed solution.printDetail
If
TRUE
(the default), information is printed. If an integeri
greater then one, information is printed at veryi
th iteration.printBar
If
TRUE
(default isFALSE
), atxtProgressBar
(from package utils) is printed. The progress bar is not shown ifprintDetail
is an integer greater than 1.storeF
if
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
Default is
FALSE
. IfTRUE
, the solutions (ie, decision variables) in every generation are stored and returned in listxlist
(see Value section below). To check, for instance, the current solution at the end of thei
th generation, retrievexlist[[c(2L, i)]]
.classify
Logical; default is
FALSE
. IfTRUE
, the result will have a class attributeSAopt
attached.OF.target
Numeric; when specified, the algorithm will stop when an objectivefunction value as low as
OF.target
(or lower) is achieved. This is useful when an optimal objectivefunction value is known: the algorithm will then stop and not waste time searching for a better solution.
At the minimum, algo
needs to contain an initial solution
x0
and a neighbour
function.
The total number of iterations equals algo$nT
times
algo$nS
(plus possibly algo$nD
).
Value
SAopt
returns a list with five components:
xbest 
the solution 
OFvalue 
objective function value of the solution, ie,

Fmat 
if 
xlist 
if 
initial.state 
the value of 
If algo$classify
was set to TRUE
, the resulting list
will have a class attribute TAopt
.
Note
If the ...
argument is used, then all the objects passed
with ...
need to go into the objective function and the
neighbourhood function. It is recommended to collect all information
in a list myList
and then write OF
and neighbour
so that they are called as OF(x, myList)
and neighbour(x,
myList)
. Note that x
need not be a vector but can be any data
structure (eg, a matrix
or a list
).
Using an initial and final temperature of zero means that
SA will be equivalent to a Local Search. The function
LSopt
may be preferred then because of smaller
overhead.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Kirkpatrick, S., Gelatt, C.D. and Vecchi, M.P. (1983). Optimization with Simulated Annealing. Science. 220 (4598), 671–680.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## Aim: given a matrix x with n rows and 2 columns,
## divide the rows of x into two subsets such that
## in one subset the columns are highly correlated,
## and in the other lowly (negatively) correlated.
## constraint: a single subset should have at least 40 rows
## create data with specified correlation
n < 100L
rho < 0.7
C < matrix(rho, 2L, 2L); diag(C) < 1
x < matrix(rnorm(n * 2L), n, 2L) %*% chol(C)
## collect data
data < list(x = x, n = n, nmin = 40L)
## a random initial solution
x0 < runif(n) > 0.5
## a neighbourhood function
neighbour < function(xc, data) {
xn < xc
p < sample.int(data$n, size = 1L)
xn[p] < abs(xn[p]  1L)
# reject infeasible solution
c1 < sum(xn) >= data$nmin
c2 < sum(xn) <= (data$n  data$nmin)
if (c1 && c2) res < xn else res < xc
as.logical(res)
}
## check (should be 1 FALSE and n1 TRUE)
x0 == neighbour(x0, data)
## objective function
OF < function(xc, data)
abs(cor(data$x[xc, ])[1L, 2L]  cor(data$x[!xc, ])[1L, 2L])
## check
OF(x0, data)
## check
OF(neighbour(x0, data), data)
## plot data
par(mfrow = c(1,3), bty = "n")
plot(data$x,
xlim = c(3,3), ylim = c(3,3),
main = "all data", col = "darkgreen")
## *Local Search*
algo < list(nS = 3000L,
neighbour = neighbour,
x0 = x0,
printBar = FALSE)
sol1 < LSopt(OF, algo = algo, data=data)
sol1$OFvalue
## *Simulated Annealing*
algo$nT < 10L
algo$nS < ceiling(algo$nS/algo$nT)
sol < SAopt(OF, algo = algo, data = data)
sol$OFvalue
c1 < cor(data$x[ sol$xbest, ])[1L, 2L]
c2 < cor(data$x[!sol$xbest, ])[1L, 2L]
lines(data$x[ sol$xbest, ], type = "p", col = "blue")
plot(data$x[ sol$xbest, ], col = "blue",
xlim = c(3, 3), ylim = c(3, 3),
main = paste("subset 1, corr.", format(c1, digits = 3)))
plot(data$x[!sol$xbest, ], col = "darkgreen",
xlim = c(3,3), ylim = c(3,3),
main = paste("subset 2, corr.", format(c2, digits = 3)))
## compare LS/SA
par(mfrow = c(1, 1), bty = "n")
plot(sol1$Fmat[ , 2L],type = "l", ylim=c(1.5, 0.5),
ylab = "OF", xlab = "Iterations")
lines(sol$Fmat[ , 2L],type = "l", col = "blue")
legend(x = "topright", legend = c("LS", "SA"),
lty = 1, lwd = 2, col = c("black", "blue"))
Download Robert Shiller's Data
Description
Download the data provided by Robert Shiller and transform them into a data frame.
Usage
Shiller(dest.dir,
url = "http://www.econ.yale.edu/~shiller/data/ie_data.xls")
Arguments
dest.dir 
character: a path to a directory 
url 
the data URL 
Details
The function downloads US stockmarket data provided by Robert Shiller which he used in his book ‘Irrational Exhuberance’. Since the data are provided in Excel format, package readxl is required.
The downloaded Excel gets a date prefix (today in
format YYYYMMDD
) and is stored in directory
dest.dir
. Before any download is attempted,
the function checks whether a file with today's
prefix exist in dest.dir
; if yes, the file is
used.
Value
a data.frame
:
Date 
end of month 
Price 
numeric 
Dividend 
numeric 
Earnings 
numeric 
CPI 
numeric 
Long Rate 
numeric 
Real Price 
numeric 
Real Dividend 
numeric 
Real Earnings 
numeric 
CAPE 
numeric 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Shiller, R.J. (2015) Irrational Exhuberance. Princeton University Press. 3rd edition.
See Also
Examples
## Not run:
archive.dir < "~/Downloads/Shiller"
if (!dir.exists(archive.dir))
dir.create(archive.dir)
Shiller(archive.dir)
## End(Not run)
Display Code Examples
Description
Display the code examples from ‘Numerical Methods and Optimization and Finance’.
Usage
showExample(file = "", chapter = NULL, showfile = TRUE,
includepaths = FALSE, edition = 2, search,
..., ignore.case = TRUE)
showChapterNames(edition = 2)
Arguments
file 
a character vector of length one. See Details. 
chapter 
optional: a character vector of length one, giving the chapter name
(see Details), or an integer, indicating a chapter number. Default
is 
showfile 
Should the file be displayed with 
includepaths 
Should the file paths be displayed? Defaults to 
... 
Arguments passed to 
edition 
an integer: 
search 
a regular expression: search in the code files. Not supported yet. 
ignore.case 
passed to 
Details
showExample
matches the specified file
argument against the available file names via
grepl(file, all.filenames, ignore.case =
ignore.case, ...)
. If chapter
is specified, a
second match is performed, grepl(chapter,
all.chapternames, ignore.case = ignore.case, ...)
. The
chapternames
are those in the book (e.g.,
‘Modeling dependencies’). The selected files
are then those for which file name and chapter name
could be matched.
Value
showExample
returns a data.frame
of at least two
character vectors, Chapter and File. If includepaths
is
TRUE
, Paths are included. If no file is found, the
data.frame
has zero rows. If a single file is identified
and showfile
is TRUE
, the function has the side effect
of displaying that file.
showChapterNames
returns a character vector: the
names of the book's chapters.
Note
The behaviour of the function changed slightly with
version 2.0 to accommodate the code examples of the
second edition of the book. Specifically, the
function gained an argument edition
, which
defaults to 2
. Also, the default for
ignore.case
was changed to TRUE
. To
get back the old behaviour of the function, set
edition
to 1
and ignore.case
to
FALSE
.
The code files can also be downloaded from https://gitlab.com/NMOF .
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2011) Numerical Methods and Optimization in Finance. Elsevier. doi:10.1016/C20090305693
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
## list all files
showExample() ## 2nd edition is default
showExample(edition = 1)
## list specific files
showExample("Appendix")
showExample("Backtesting")
showExample("Heuristics")
showExample("tutorial") ## matches against filename
showExample(chapter = 13)
showExample(chapter = "tutorial")
## show where a file is installed
showExample(chapter = "portfolio", includepaths = TRUE)
## first edition
showExample("equations.R", edition = 1)
showExample("example", chapter = "portfolio", edition = 1)
showExample("example", chapter = 13, edition = 1)
showExample("example", chapter = showChapterNames(1)[13L], edition = 1)
ThresholdAccepting Information
Description
The function can be called from the objective and neighbourhood
function during a run of TAopt
; it provides information
such as the current iteration, the current solution, etc.
Usage
TA.info(n = 0L)
Arguments
n 
generational offset; see Details. 
Details
This function is still experimental.
The function can be called in the neighbourhood function or the
objective function during a run of TAopt
. It evaluates
to a list with the state of the optimisation run, such as the current
iteration.
TA.info
relies on parent.frame
to retrieve its
information. If the function is called within another function in the
neighbourhood or objective function, the argument n
needs to be
increased.
Value
A list
OF.sampling 
logical: if 
iteration 
current iteration 
step 
current step (i.e. for a given threshold) 
threshold 
current threshold (the number, not the value) 
xbest 
the best solution found so far 
OF.xbest 
objective function value of best solution 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
### MINIMAL EXAMPLE for TAopt
## objective function evaluates to a constant
fun < function(x)
0
## neighbourhood function does not even change the solution,
## but it reports information
nb < function(x) {
tmp < TA.info()
cat("current threshold ", tmp$threshold,
" current step ", tmp$step,
" current iteration ", tmp$iteration, "\n")
x
}
## run TA
algo < list(nS = 5,
nT = 2,
nD = 3,
x0 = rep(0, 5),
neighbour = nb,
printBar = FALSE,
printDetail = FALSE)
ignore < TAopt(fun, algo)
## printed output:
## current threshold NA  current step 1  current iteration NA
## current threshold NA  current step 2  current iteration NA
## current threshold NA  current step 3  current iteration NA
## current threshold 1  current step 1  current iteration 1
## current threshold 1  current step 2  current iteration 2
## current threshold 1  current step 3  current iteration 3
## current threshold 1  current step 4  current iteration 4
## current threshold 1  current step 5  current iteration 5
## current threshold 2  current step 1  current iteration 6
## current threshold 2  current step 2  current iteration 7
## current threshold 2  current step 3  current iteration 8
## current threshold 2  current step 4  current iteration 9
## current threshold 2  current step 5  current iteration 10
Optimisation with Threshold Accepting
Description
The function implements the Threshold Accepting algorithm.
Usage
TAopt(OF, algo = list(), ...)
Arguments
OF 
The objective function, to be minimised. Its first argument
needs to be a solution 
algo 
A list of settings for the algorithm. See Details. 
... 
other variables passed to 
Details
Threshold Accepting (TA) changes an initial solution
iteratively; the algorithm stops after a fixed number of
iterations. Conceptually, TA consists of a loop than runs
for a number of iterations. In each iteration, a current solution
xc
is changed through a function algo$neighbour
. If this
new (or neighbour) solution xn
is not worse than xc
, ie,
if OF(xn,...) <= OF(xc,...)
, then xn
replaces
xc
. If xn
is worse, it still replaces xc
as long
as the difference in ‘quality’ between the two solutions is
less than a threshold tau
; more precisely, as long as
OF(xn,...)  tau <= OF(xc,...)
. Thus, we also accept a new
solution that is worse than its predecessor; just not too much
worse. The threshold is typically decreased over the course of the
optimisation. For zero thresholds TA becomes a stochastic local
search.
The thresholds can be passed through the list algo
(see
below). Otherwise, they are automatically computed through the
procedure described in Gilli et al. (2006). When the thresholds are
created automatically, the final threshold is always zero.
The list algo
contains the following items.
nS
The number of steps per threshold. The default is 1000; but this setting depends very much on the problem.
nT
The number of thresholds. Default is 10; ignored if
algo$vT
is specified.nI

Total number of iterations, with default
NULL
. If specified, it will overridenS
withceiling(nI/nT)
. Using this option makes it easier to compare and switch between functionsLSopt
,TAopt
andSAopt
. nD
The number of random steps to compute the threshold sequence. Defaults to 2000. Only used if
algo$vT
isNULL
.q
The highest quantile for the threshold sequence. Defaults to 0.5. Only used if
algo$vT
isNULL
. Ifq
is zero,TAopt
will run withalgo$nT
zerothresholds (ie, like a Local Search).x0
The initial solution. If this is a function, it will be called once without arguments to compute an initial solution, ie,
x0 < algo$x0()
. This can be useful when the routine is called in a loop of restarts, and each restart is to have its own starting value.vT
The thresholds. A numeric vector. If
NULL
(the default),TAopt
will computealgo$nT
thresholds. Passing threshold can be useful when similar problems are handled. Then the time to sample the objective function to compute the thresholds can be saved (ie, we savealgo$nD
function evaluations). If the thresholds are computed andalgo$printDetail
isTRUE
, the time required to evaluate the objective function will be measured and an estimate for the remaining computing time is issued. This estimate is often very crude.neighbour
The neighbourhood function, called as
neighbour(x, ...)
. Its first argument must be a solutionx
; it must return a changed solution.printDetail
If
TRUE
(the default), information is printed. If an integeri
greater then one, information is printed at veryi
th iteration.printBar
If
TRUE
(default isFALSE
), atxtProgressBar
(from package utils) is printed. The progress bar is not shown ifprintDetail
is an integer greater than 1.scale
The thresholds are multiplied by
scale
. Default is 1.drop0
When thresholds are computed, should zero values be dropped from the sample of objectivefunction values? Default is
FALSE
.stepUp
Defaults to
0
. If an integer greater than zero, then the thresholds are recycled, ie,vT
is replaced byrep(vT, algo$stepUp + 1)
(and the number of thresholds will be increased byalgo$nT
timesalgo$stepUp
). This option works for supplied as well as computed thresholds. Practically, this will have the same effect as restarting from a returned solution. (In Simulated Annealing, this strategy goes by the name of ‘reheating’.)thresholds.only
Defaults to
FALSE
. IfTRUE
, compute only threshold sequence, but do not actually run TA.storeF
if
TRUE
(the default), the objective function values for every solution in every generation are stored and returned as matrixFmat
.storeSolutions
Default is
FALSE
. IfTRUE
, the solutions (ie, decision variables) in every generation are stored and returned in listxlist
(see Value section below). To check, for instance, the current solution at the end of thei
th generation, retrievexlist[[c(2L, i)]]
.classify
Logical; default is
FALSE
. IfTRUE
, the result will have a class attributeTAopt
attached. This feature is experimental: the supported methods (plot, summary) may change without warning.OF.target
Numeric; when specified, the algorithm will stop when an objectivefunction value as low as
OF.target
(or lower) is achieved. This is useful when an optimal objectivefunction value is known: the algorithm will then stop and not waste time searching for a better solution.
At the minimum, algo
needs to contain an initial solution
x0
and a neighbour
function.
The total number of iterations equals algo$nT
times
(algo$stepUp + 1)
times algo$nS
(plus possibly
algo$nD
).
Value
TAopt
returns a list with four components:
xbest 
the solution 
OFvalue 
objective function value of the solution, ie,

Fmat 
if 
xlist 
if 
initial.state 
the value of 
If algo$classify
was set to TRUE
, the resulting list
will have a class attribute TAopt
.
Note
If the ...
argument is used, then all the objects passed
with ...
need to go into the objective function and the
neighbourhood function. It is recommended to collect all information
in a list myList
and then write OF
and neighbour
so that they are called as OF(x, myList)
and neighbour(x,
myList)
. Note that x
need not be a vector but can be any data
structure (eg, a matrix
or a list
).
Using thresholds of size 0 makes TA run as a Local Search. The
function LSopt
may be preferred then because of smaller
overhead.
Author(s)
Enrico Schumann
References
Dueck, G. and Scheuer, T. (1990) Threshold Accepting. A General Purpose Optimization Algorithm Superior to Simulated Annealing. Journal of Computational Physics. 90 (1), 161–175.
Dueck, G. and Winker, P. (1992) New Concepts and Algorithms for Portfolio Choice. Applied Stochastic Models and Data Analysis. 8 (3), 159–178.
Gilli, M., Këllezi, E. and Hysi, H. (2006) A DataDriven Optimization Heuristic for Downside Risk Minimization. Journal of Risk. 8 (3), 1–18.
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Moscato, P. and Fontanari, J.F. (1990). Stochastic Versus Deterministic Update in Simulated Annealing. Physics Letters A. 146 (4), 204–208.
Schumann, E. (2012) Remarks on 'A comparison of some heuristic optimization methods'. http://enricoschumann.net/R/remarks.htm
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Winker, P. (2001). Optimization Heuristics in Econometrics: Applications of Threshold Accepting. Wiley.
See Also
LSopt
, restartOpt
. Simulated Annealing
is implemented in function SAopt
.
Package neighbours (also on CRAN) offers helpers
for creating neighbourhood functions.
Examples
## Aim: given a matrix x with n rows and 2 columns,
## divide the rows of x into two subsets such that
## in one subset the columns are highly correlated,
## and in the other lowly (negatively) correlated.
## constraint: a single subset should have at least 40 rows
## create data with specified correlation
n < 100L
rho < 0.7
C < matrix(rho, 2L, 2L); diag(C) < 1
x < matrix(rnorm(n * 2L), n, 2L) %*% chol(C)
## collect data
data < list(x = x, n = n, nmin = 40L)
## a random initial solution
x0 < runif(n) > 0.5
## a neighbourhood function
neighbour < function(xc, data) {
xn < xc
p < sample.int(data$n, size = 1L)
xn[p] < abs(xn[p]  1L)
# reject infeasible solution
c1 < sum(xn) >= data$nmin
c2 < sum(xn) <= (data$n  data$nmin)
if (c1 && c2) res < xn else res < xc
as.logical(res)
}
## check (should be 1 FALSE and n1 TRUE)
x0 == neighbour(x0, data)
## objective function
OF < function(xc, data)
abs(cor(data$x[xc, ])[1L, 2L]  cor(data$x[!xc, ])[1L, 2L])
## check
OF(x0, data)
## check
OF(neighbour(x0, data), data)
## plot data
par(mfrow = c(1,3), bty = "n")
plot(data$x,
xlim = c(3,3), ylim = c(3,3),
main = "all data", col = "darkgreen")
## *Local Search*
algo < list(nS = 3000L,
neighbour = neighbour,
x0 = x0,
printBar = FALSE)
sol1 < LSopt(OF, algo = algo, data=data)
sol1$OFvalue
## *Threshold Accepting*
algo$nT < 10L
algo$nS < ceiling(algo$nS/algo$nT)
sol < TAopt(OF, algo = algo, data = data)
sol$OFvalue
c1 < cor(data$x[ sol$xbest, ])[1L, 2L]
c2 < cor(data$x[!sol$xbest, ])[1L, 2L]
lines(data$x[ sol$xbest, ], type = "p", col = "blue")
plot(data$x[ sol$xbest, ], col = "blue",
xlim = c(3,3), ylim = c(3,3),
main = paste("subset 1, corr.", format(c1, digits = 3)))
plot(data$x[!sol$xbest, ], col = "darkgreen",
xlim = c(3,3), ylim = c(3,3),
main = paste("subset 2, corr.", format(c2, digits = 3)))
## compare LS/TA
par(mfrow = c(1,1), bty = "n")
plot(sol1$Fmat[ ,2L],type="l", ylim=c(1.5,0.5),
ylab = "OF", xlab = "iterations")
lines(sol$Fmat[ ,2L],type = "l", col = "blue")
legend(x = "topright",legend = c("LS", "TA"),
lty = 1, lwd = 2,col = c("black", "blue"))
Classical Test Functions for Unconstrained Optimisation
Description
A number of functions that have been suggested in the literature as benchmarks for unconstrained optimisation.
Usage
tfAckley(x)
tfEggholder(x)
tfGriewank(x)
tfRastrigin(x)
tfRosenbrock(x)
tfSchwefel(x)
tfTrefethen(x)
Arguments
x 
a numeric vector of arguments. See Details. 
Details
All functions take as argument only one variable, a
numeric vector x
whose length determines the
dimensionality of the problem.
The Ackley function is implemented as
\exp(1) + 20 20 \exp{\left(0.2 \sqrt{\frac{1}{n}\sum_{i=1}^n x_i^2}\right)}  \exp{\left(\frac{1}{n}\sum_{i=1}^n \cos(2 \pi x_i)\right)}\,.
The minimum function value is zero; reached at x=0
.
The Eggholder takes a twodimensional x
,
here written as x
and y
. It is defined as
(y + 47) \sin\left(\sqrt{y + \frac{x}{2} + 47}\right) 
x \sin\left(\sqrt{x  (y + 47)}\right)\,.
The minimum function value is 959.6407; reached at c(512, 404.2319)
.
The Griewank function is given by
1+\frac{1}{4000} \sum^n_{i=1} x_i^2  \prod_{i=1}^n \cos \left(\frac{x_i}{\sqrt{i}}\right)\,.
The function is minimised at x=0
; its minimum value is zero.
The Rastrigin function:
10n + \sum_{i=1}^n \left(x_i^2 10\cos(2\pi x_i)\right)\,.
The minimum function value is zero; reached at x=0
.
The Rosenbrock (or banana) function:
\sum_{i=1}^{n1}\left(100 (x_{i+1}x_i^2)^2 + (1x_i)^2\right)\,.
The minimum function value is zero; reached at x=1
.
The Schwefel function:
\sum_{i=1}^n \left(x_i \sin\left(\sqrt{x_i}\right)\right)\,.
The minimum function value (to about 8 digits) is 418.9829n
; reached at x = 420.9687
.
Trefethen's function takes a twodimensional x
(here written as x
and y
); it is defined as
\exp(\sin(50x)) + \sin(60 e^y) + \sin(70 \sin(x)) + \sin(\sin(80y))  \sin(10(x+y)) + \frac{1}{4}(x^2+y^2)\,.
The minimum function value is 3.3069; reached at c(0.0244, 0.2106)
.
Value
The objective function evaluated at x
(a numeric vector of
length one).
Warning
These test functions represent artificial problems. It is practically not too helpful to finetune a method on such functions. (That would be like memorising all the answers to a particular multiplechoice test.) The functions' main purpose is checking the numerical implementation of algorithms.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## persp for twodimensional x
## Ackley
n < 100L; surf < matrix(NA, n, n)
x1 < seq(from = 2, to = 2, length.out = n)
for (i in 1:n)
for (j in 1:n)
surf[i, j] < tfAckley(c(x1[i], x1[j]))
persp(x1, x1, surf, phi = 30, theta = 30, expand = 0.5,
col = "goldenrod1", shade = 0.2, ticktype = "detailed",
xlab = "x1", ylab = "x2", zlab = "f", main = "Ackley (f)",
border = NA)
## Trefethen
n < 100L; surf < matrix(NA, n, n)
x1 < seq(from = 10, to = 10, length.out = n)
for (i in 1:n)
for (j in 1:n)
surf[i, j] < tfTrefethen(c(x1[i], x1[j]))
persp(x1, x1, surf, phi = 30, theta = 30, expand = 0.5,
col = "goldenrod1", shade = 0.2, ticktype = "detailed",
xlab = "x1", ylab = "x2", zlab = "f", main = "Trefethen (f)",
border = NA)
Compute a Tracking Portfolio
Description
Computes a portfolio similar to a benchmark, e.g. for tracking the benchmark's performance or identifying factors.
Usage
trackingPortfolio(var, wmin = 0, wmax = 1,
method = "qp", objective = "variance", R,
ls.algo = list())
Arguments
var 
the covariance matrix: a numeric (real), symmetric matrix. The first asset is the benchmark. 
R 
a matrix of returns: each colums holds the returns of one asset; each rows holds the returns for one observation. The first asset is the benchmark. 
wmin 
numeric: a lower bound on weights. May also be a vector that holds specific bounds for each asset. 
wmax 
numeric: an upper bound on weights. May also be a vector that holds specific bounds for each asset. 
method 
character. Currently, 
objective 
character. Currently, 
ls.algo 
a list of named elements, for settings for
method ‘ 
Details
With method "qp"
, the function uses
solve.QP
from package
quadprog. Because of the algorithm that
solve.QP
uses, var
has to
be positive definite (i.e. must be of full rank).
With method "ls"
, the function uses
LSopt
. Settings can be passed via
ls.algo
, which corresponds to
LSopt
's argument algo
. Default
settings are 2000 iterations and printBar
,
printDetail
set to FALSE
.
R
is needed only when objective
is
"sum.of.squares"
or method
is
‘ls
’. (See Examples.)
Value
a numeric vector (the portfolio weights)
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Schumann, E. (2020) Returnbased tracking portfolios. http://enricoschumann.net/notes/returnbasedtrackingportfolios.html
Sharpe, W. F. (1992) Asset Allocation: Management Style and Performance Measurement. Journal of Portfolio Management. 18 (2), 7–19. https://web.stanford.edu/~wfsharpe/art/sa/sa.htm
See Also
Examples
if (requireNamespace("quadprog")) {
ns < 120
R < randomReturns(na = 1 + 20,
ns = ns,
sd = 0.03,
mean = 0.005,
rho = 0.7)
var < cov(R)
sol.qp < trackingPortfolio(var, wmax = 0.4)
sol.ls < trackingPortfolio(var = var, R = R, wmax = 0.4, method = "ls")
data.frame(QP = round(100*sol.qp, 1),
LS = round(100*sol.ls, 1))
sol.qp < trackingPortfolio(var, R = R, wmax = 0.4,
objective = "sum.of.squares")
sol.ls < trackingPortfolio(var = var, R = R, wmax = 0.4, method = "ls",
objective = "sum.of.squares")
data.frame(QP = round(100*sol.qp, 1),
LS = round(100*sol.ls, 1))
## same as 'sol.qp' above
sol.qp.R < trackingPortfolio(R = R,
wmax = 0.4,
objective = "sum.of.squares")
sol.qp.var < trackingPortfolio(var = crossprod(R),
wmax = 0.4,
objective = "variance")
## ==> should be the same
all.equal(sol.qp.R, sol.qp.var)
}
Pricing PlainVanilla Bonds
Description
Calculate the theoretical price and yieldtomaturity of a list of cashflows.
Usage
vanillaBond(cf, times, df, yields)
ytm(cf, times, y0 = 0.05, tol = 1e05, maxit = 1000L, offset = 0)
duration(cf, times, yield, modified = TRUE, raw = FALSE)
convexity(cf, times, yield, raw = FALSE)
Arguments
cf 
Cashflows; a numeric vector or a matrix. If a matrix, cashflows should be arranged in rows; timestopayment correspond to columns. 
times 
timestopayment; a numeric vector 
df 
discount factors; a numeric vector 
yields 
optional (instead of discount factors); zero yields to compute discount factor; if of length one, a flat zero curve is assumed 
yield 
numeric vector of length one (both duration and convexity assume a flat yield curve) 
y0 
starting value 
tol 
tolerance 
maxit 
maximum number of iterations 
offset 
numeric: a ‘base’ rate over which to compute the yield to maturity. See Details and Examples. 
modified 
logical: return modified duration? (default 
raw 
logical: default 
Details
vanillaBond
computes the present value of a vector of
cashflows; it may thus be used to evaluate not just bonds but any
instrument that can be reduced to a deterministic set of cashflows.
ytm
uses Newton's method to compute the yieldtomaturity of a
bond (a.k.a. internal interest rate). When used with a bond, the initial
outlay (i.e. the bonds dirty price) needs be included in the vector of
cashflows. For a coupon bond, a good starting value y0
is
the coupon divided by the dirty price of the bond.
An offset
can be specified either as a single number or as a
vector of zero rates. See Examples.
Value
numeric
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
See Also
Examples
## ytm
cf < c(5, 5, 5, 5, 5, 105) ## cashflows
times < 1:6 ## maturities
y < 0.0127 ## the "true" yield
b0 < vanillaBond(cf, times, yields = y)
cf < c(b0, cf); times < c(0, times)
ytm(cf, times)
## ... with offset
cf < c(5, 5, 5, 5, 5, 105) ## cashflows
times < 1:6 ## maturities
y < 0.02 + 0.01 ## riskfree 2% + riskpremium 1%
b0 < vanillaBond(cf, times, yields = y)
cf < c(b0, cf); times < c(0, times)
ytm(cf, times, offset = 0.02) ## ... only the riskpremium
cf < c(5, 5, 5, 5, 5, 105) ## cashflows
times < 1:6 ## maturities
y < NS(c(6,9,10,5)/100, times) ## riskpremium 1%
b0 < vanillaBond(cf, times, yields = y + 0.01)
cf < c(b0, cf); times < c(0, times)
ytm(cf, times, offset = c(0,y)) ## ... only the riskpremium
## bonds
cf < c(5, 5, 5, 5, 5, 105) ## cashflows
times < 1:6 ## maturities
df < 1/(1+y)^times ## discount factors
all.equal(vanillaBond(cf, times, df),
vanillaBond(cf, times, yields = y))
## ... using NelsonSiegel
vanillaBond(cf, times, yields = NS(c(0.03,0,0,1), times))
## several bonds
## cashflows are numeric vectors in a list 'cf',
## timestopayment are are numeric vectors in a
## list 'times'
times < list(1:3,
1:4,
0.5 + 0:5)
cf < list(c(6, 6, 106),
c(4, 4, 4, 104),
c(2, 2, 2, 2, 2, 102))
alltimes < sort(unique(unlist(times)))
M < array(0, dim = c(length(cf), length(alltimes)))
for (i in seq_along(times))
M[i, match(times[[i]], alltimes)] < cf[[i]]
rownames(M) < paste("bond.", 1:3, sep = "")
colnames(M) < format(alltimes, nsmall = 1)
vanillaBond(cf = M, times = alltimes, yields = 0.02)
## duration/convexity
cf < c(5, 5, 5, 5, 5, 105) ## cashflows
times < 1:6 ## maturities
y < 0.0527 ## yield to maturity
d < 0.001 ## change in yield (+10 bp)
vanillaBond(cf, times, yields = y + d)  vanillaBond(cf, times, yields = y)
duration(cf, times, yield = y, raw = TRUE) * d
duration(cf, times, yield = y, raw = TRUE) * d +
convexity(cf, times, yield = y, raw = TRUE)/2 * d^2
Pricing PlainVanilla (European and American) and Barrier Options (European)
Description
Functions to calculate the theoretical prices and (some) Greeks for plainvanilla and barrier options.
Usage
vanillaOptionEuropean(S, X, tau, r, q, v, tauD = 0, D = 0,
type = "call", greeks = TRUE,
model = NULL, ...)
vanillaOptionAmerican(S, X, tau, r, q, v, tauD = 0, D = 0,
type = "call", greeks = TRUE, M = 101)
vanillaOptionImpliedVol(exercise = "european", price,
S, X, tau, r, q = 0,
tauD = 0, D = 0,
type = "call",
M = 101,
uniroot.control = list(),
uniroot.info = FALSE)
barrierOptionEuropean(S, X, H, tau, r, q = 0, v, tauD = 0, D = 0,
type = "call",
barrier.type = "downin",
rebate = 0,
greeks = FALSE,
model = NULL, ...)
Arguments
S 
spot 
X 
strike 
H 
barrier 
tau 
timetomaturity in years 
r 
riskfree rate 
q 
continuous dividend yield, see Details. 
v 
variance (volatility squared) 
tauD 
vector of timestodividends in years. Only dividends with

D 
vector of dividends (in currency units); default is no dividends. 
type 
call or put; default is call. 
barrier.type 
string: combination of 
rebate 
currently not implemented 
greeks 
compute Greeks? Defaults to 
model 
what model to use to value the option. Default is 
... 
parameters passed to pricing model 
M 
number of time steps in the tree 
exercise 

price 
numeric; the observed price to be recovered through choice of volatility. 
uniroot.control 
A list. If there are elements named

uniroot.info 
logical; default is 
Details
For European options the formula of Messrs Black, Scholes and
Merton is used. It can be used for equities (set q
equal
to the dividend yield), futures (Black, 1976; set q
equal
to r
), currencies (Garman and Kohlhagen, 1983; set
q
equal to the foreign riskfree rate). For futurestyle
options (e.g. options on the German Bund future), set q
and r
equal to zero.
The Greeks are provided in their raw (‘textbook’) form with only one exception: Theta is made negative. For practical use, the other Greeks are also typically adjusted: Theta is often divided by 365 (or some other yearly day count); Vega and Rho are divided by 100 to give the sensitivity for one percentagepoint move in volatility/the interest rate. Raw Gamma is not much use if not adjusted for the actual move in the underlier.
For European options the Greeks are computed through the respective analytic expressions. For American options only Delta, Gamma and Theta are computed because they can be directly obtained from the binomial tree; other Greeks need to be computed through a finite difference (see Examples).
For the Europeantype options, the function understands vectors of inputs, except for dividends. American options are priced via a CoxRossRubinstein tree; no vectorisation is implemented here.
The implied volatility is computed with uniroot
from the stats package (the default search interval is
c(0.00001, 2)
; it can be changed through
uniroot.control
).
Dividends (D
) are modelled via the escroweddividend
model.
Value
Returns the price (a numeric vector of length one) if greeks
is
FALSE
, else returns a list.
Note
If greeks
is TRUE
, the function will return a list with
named elements (value
, delta
and so on). Prior to version
0.263, the first element of this list was named price
.
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Haug, E. (2007) The Complete Guide to Option Pricing Formulas. 2nd edition. McGrawHill.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
S < 100; X < 100; tau < 1; r < 0.02; q < 0.06; vol < 0.3
unlist(vanillaOptionEuropean(S, X, tau, r, q, vol^2, type = "put"))
S < 100; X < 110; tau < 1; r < 0.1; q < 0.06; vol < 0.3; type < "put"
unlist(vanillaOptionAmerican(S, X, tau, r, q, vol^2, type = type,
greeks = TRUE))
## compute rho for 1% move
h < 0.01
(vanillaOptionAmerican(S, X, tau, r + h, q, vol^2,
type = type, greeks = FALSE) 
vanillaOptionAmerican(S, X, tau, r, q, vol^2,
type = type, greeks = FALSE)) / (h*100)
## compute vega for 1% move
h < 0.01
(vanillaOptionAmerican(S, X, tau, r, q,(vol + h)^2,
type = type, greeks = FALSE) 
vanillaOptionAmerican(S, X, tau, r, q, vol^2,
type = type, greeks = FALSE)) / (h*100)
S < 100; X < 100
tau < 1; r < 0.05; q < 0.00
D < c(1,2); tauD < c(0.3,.6)
type < "put"
v < 0.245^2 ## variance, not volatility
p < vanillaOptionEuropean(S = S, X = X, tau, r, q, v = v,
tauD = tauD, D = D, type = type, greeks = FALSE)
vanillaOptionImpliedVol(exercise = "european", price = p,
S = S, X = X, tau = tau, r = r, q = q, tauD = tauD, D = D, type = type)
p < vanillaOptionAmerican(S = S, X = X, tau, r, q, v = v,
tauD = tauD, D = D, type = type, greeks = FALSE)
vanillaOptionImpliedVol(exercise = "american", price = p,
S = S, X = X, tau = tau, r = r, q = q, tauD = tauD, D = D, type =
type, uniroot.control = list(interval = c(0.01, 0.5)))
## compute implied q
S < 100; X < 100
tau < 1; r < 0.05; q < 0.072
v < 0.22^2 ## variance, not volatility
call < vanillaOptionEuropean(S=S, X = X, tau=tau, r=r, q=q, v=v,
type = "call", greeks = FALSE)
put < vanillaOptionEuropean(S=S, X = X, tau=tau, r=r, q=q, v=v,
type = "put", greeks = FALSE)
# ... the simple way
(log(call + X * exp(tau*r)  put)  log(S)) / tau
# ... the complicated way :)
volDiffCreate < function(exercise, call, put, S, X, tau, r) {
f < function(q) {
cc < vanillaOptionImpliedVol(exercise = exercise, price = call,
S = S, X = X, tau = tau, r = r, q = q, type = "call")
pp < vanillaOptionImpliedVol(exercise = exercise, price = put,
S = S, X = X, tau = tau, r = r, q = q, type = "put")
abs(cc  pp)
}
f
}
f < volDiffCreate(exercise = "european",
call = call, put = put, S = S, X = X, tau = tau, r)
optimise(f,interval = c(0, 0.2))$minimum
##
S < 100; X < 100
tau < 1; r < 0.05; q < 0.072
v < 0.22^2 ## variance, not volatility
vol < 0.22
vanillaOptionEuropean(S=S, X = X, tau=tau, r=r, q=q, v=v, ## with variance
type = "call", greeks = FALSE)
vanillaOptionEuropean(S=S, X = X, tau=tau, r=r, q=q, vol=vol, ## with vol
type = "call", greeks = FALSE)
vanillaOptionEuropean(S=S, X = X, tau=tau, r=r, q=q, vol=vol, ## with vol
type = "call", greeks = FALSE, v = 0.2^2)
Contract Value of Australian Government Bond Future
Description
Compute the contract value of an Australian governmentbond future from its quoted price.
Usage
xtContractValue(quoted.price, coupon, do.round = TRUE)
xtTickValue(quoted.price, coupon, do.round = TRUE)
Arguments
quoted.price 
The price, as in 
coupon 
numeric; should be 6, not 0.06 
do.round 
If 
Details
Australian governmentbond futures, traded at the
Australian Securities Exchange (asx), are
quoted as 100  yield
. The function computes
the actual contract value from the quoted price.
xtTickValue
computes the tick value via a
central difference.
Value
A numeric vector.
Author(s)
Enrico Schumann
References
https://www.rba.gov.au/mktoperations/resources/technotes/pricingformulae.html
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
Examples
quoted.price < 99
coupon < 6
xtContractValue(quoted.price, coupon)
xtTickValue(quoted.price, coupon)
## convexity
quoted.price < seq(90, 100, by = 0.1)
plot(100  quoted.price,
xtContractValue(quoted.price, coupon),
xlab = "Yield", ylab = "Contract value")
Integration of Gausstype
Description
Compute nodes and weights for Gauss integration.
Usage
xwGauss(n, method = "legendre")
changeInterval(nodes, weights, oldmin, oldmax, newmin, newmax)
Arguments
n 
number of nodes 
method 
character. default is 
nodes 
the nodes (a numeric vector) 
weights 
the weights (a numeric vector) 
oldmin 
the minimum of the interval (typically as tabulated) 
oldmax 
the maximum of the interval (typically as tabulated) 
newmin 
the desired minimum of the interval 
newmax 
the desired maximum of the interval 
Details
xwGauss
computes nodes and weights for integration for the
interval 1 to 1. It uses the method of Golub and Welsch (1969).
changeInterval
is a utility that transforms nodes and weights
to an arbitrary interval.
Value
a list with two elements
weights 
a numeric vector 
nodes 
a numeric vector 
Author(s)
Enrico Schumann
References
Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017001621X
Golub, G.H. and Welsch, J.H. (1969). Calculation of Gauss Quadrature Rules. Mathematics of Computation, 23(106), pp. 221–230+s1–s10.
Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). http://enricoschumann.net/NMOF.htm#NMOFmanual
See Also
Examples
## examples from Gilli/Maringer/Schumann (2019), ch. 17
## a test function
f1 < function(x) exp(x)
m < 5; a < 0; b < 5
h < (b  a)/m
## rectangular rule  left
w < h; k < 0:(m1); x < a + k * h
sum(w * f1(x))
## rectangular rule  right
w < h; k < 1:m ; x < a + k * h
sum(w * f1(x))
## midpoint rule
w < h; k < 0:(m1); x < a + (k + 0.5)*h
sum(w * f1(x))
## trapezoidal rule
w < h
k < 1:(m1)
x < c(a, a + k*h, b)
aux < w * f1(x)
sum(aux)  (aux[1] + aux[length(aux)])/2
## R's integrate (from package stats)
integrate(f1, lower = a,upper = b)
## GaussLegendre
temp < xwGauss(m)
temp < changeInterval(temp$nodes, temp$weights,
oldmin = 1, oldmax = 1, newmin = a, newmax = b)
x < temp$nodes; w < temp$weights
sum(w * f1(x))