Package 'NMOF'

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 general-purpose 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 mean-variance 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 major-minor-patch. 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 https://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

See resampleC and mc.

Data

See bundData, fundData and optionData.

Author(s)

Enrico Schumann

Maintainer: Enrico Schumann <[email protected]>

References

Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. doi:10.1016/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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.0-0", package = "NMOF")

## vignettes
vignette(package = "NMOF")
nss <- vignette("DEnss", package = "NMOF")
print(nss)
edit(nss)

## _book_ websites
browseURL("https://nmof.net")
browseURL("https://enricoschumann.net/NMOF/")

## _package_ websites
browseURL("https://enricoschumann.net/R/packages/NMOF/")
browseURL("https://cran.r-project.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(".*\\(([0-9]+) checks?\\).*", "\\1",
           test.rep[grep("\\(\\d+ checks?\\)", test.rep)])
message("Number of unit tests: ", sum(as.numeric(nt)))

---

Approximate Total Return of Bond

Description

Approximate the total return of a bond by its current yield, duration and convexity.

Usage

approxBondReturn(yield, tm, n = 2, scale = 1/250, pad = NULL)

Arguments

yield	a numeric vector
tm	a numeric vector: time-to-maturity
n	number of coupon payments per period
scale	how to scale yield; see Details
pad	how to pad the first observation: NULL (default) means to drop it; useful alternatives are NA or 0

Details

The function approximates the total return of a bond investor, based on changes in yield. The computation is based on a Taylor-series expansion. See the references, in particular concerning the shortcomings of the approximation:

1. approximation is based on par yield

2. it relies on yield alone, so does not take into account defaults; so for indices, the approximation should only used for issuers without defaults

Value

a numeric vector, with attributes duration and convexity

Note

Package treasuryTR implements the method as well.

Author(s)

Enrico Schumann

References

Swinkels, L. (2019). Treasury Bond Return Data Starting in 1962. Data. 4 (3).

Tuckman, B. and Serrat, A. (2012). Fixed Income Securities – Tools for Today's Markets. 3rd edition. Wiley.

Examples

yield0 <- 0.05
tm <- 20
cf <- c(rep(5, tm-1), 105)
duration(cf, 1:tm, yield0)

approxBondReturn(yield = c(yield0, 0.05), tm = tm, n = 1)
## ==> no price change, current yield is earned

approxBondReturn(yield = c(yield0, 0.04), tm = tm, n = 1)
## ==> current yield + price changed is earned

---

Zero-Bracketing

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 fun(x, ...) with x being a numeric vector
interval	a numeric vector, containing the end-points of the interval to be searched
...	further arguments passed to fun
lower	lower end-point. Ignored if interval is specified.
upper	upper end-point. Ignored if interval is specified.
n	the number of function evaluations. Must be at least 2 (in which case fun is evaluated only at the end-points); defaults to 20.
method	can be loop (the default), vectorised, multicore or snow. See Details.
mc.control	a list containing settings that will be passed to mclapply if method is multicore. Must be a list of named elements. See the documentation of mclapply in package parallel.
cl	default is NULL. If method is snow, this must be a cluster object or an integer (the number of cores to be used). See the documentation of packages parallel and snow.

Details

bracketing evaluates fun at equal-spaced 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 2010-05-31.

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 YYYY-MM-DD. 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.deutsche-finanzagentur.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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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	Date or character in format YYYY-MM-DD
expiry.date	Date or character in format YYYY-MM-DD
last.coupon.date	Date or character in format YYYY-MM-DD
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 cheapest-to-deliver 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Examples

## Bund-Future 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 = "2017-8-21",
           expiry.date = "2017-09-07",       ## Bund expiry
           last.coupon.date = "2017-08-15",  ## last co
           r = -0.0037,
           cf = 0.594455)   ## conversion factor (from Eurex website)

bundFutureImpliedRate(future = 164.44,
                      clean = 97.769,
                      coupon = 0,
                      trade.date = "2017-8-21",
                      expiry.date = "2017-09-07",
                      last.coupon.date = "2017-08-15",
                      cf = 0.594455)

---

Price a Plain-Vanilla Call with the Characteristic Function

Description

Price a European plain-vanilla 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 interval, tol or maxiter, these are passed to uniroot. Any other elements of the list are ignored.
uniroot.info	logical; default is FALSE. If TRUE, the function will return the information returned by uniroot. See paragraph Value below.
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.26-3, 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/C2017-0-01621-X

Heston, S.L. (1993) A Closed-Form 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). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

callHestoncf

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)

## Black-Scholes-Merton
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	risk-free rate
q	dividend rate
v0	current variance
vT	long-run variance (theta in Heston's paper)
rho	correlation between spot and variance
k	speed of mean-reversion (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 FALSE.
...	named arguments, passed to integrate
uniroot.control	A list. If there are elements named interval, tol or maxiter, these are passed to uniroot. Other elements of the list are ignored.
uniroot.info	logical; default is FALSE. If TRUE, the function will return the information returned by uniroot. See section Value below.

Details

The function computes the value of a plain vanilla European call under the Heston model. Put values can be computed through put–call-parity.

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.26-3, 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/C2017-0-01621-X

Heston, S.L. (1993) A Closed-Form 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). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

callCF, EuropeanCall

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	risk-free rate
q	dividend rate
v	variance
lambda	jump intensity
muJ	mean jump-size
vJ	variance of log jump-size
N	The number of jumps. See Details.
implVol	compute equivalent Black–Scholes–Merton volatility? Default is FALSE.

Details

The function computes the value of a plain-vanilla European call under Merton's jump–diffusion model. Put values can be computed through put–call-parity (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/C2017-0-01621-X

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). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

callCF, EuropeanCall

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)

---

Full-rank Column Subset

Description

Select a full-rank 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

repairMatrix

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)

---

Constant-Proportion Portfolio Insurance

Description

Simulate constant-proportion 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 population-based 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 and max 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 if algo$minmaxConstr is TRUE will the algorithm repair solutions outside the min and max range.

minmaxConstr

if TRUE, algo$min and algo$max are considered constraints. Default is FALSE.

pen

a penalty function. Default is NULL (no penalty).

initP

optional: the initial population. A matrix of size length(algo$min) times algo$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 to TRUE.

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 integer i greater then one, information is printed at very ith generation.

printBar

If TRUE (the default), a txtProgressBar is printed.

storeF

if TRUE (the default), the objective function values for every solution in every generation are stored and returned as matrix Fmat.

storeSolutions

default is FALSE. If TRUE, the solutions (ie, decision variables) in every generation are stored and returned as a list P in list xlist (see Value section below). To check, for instance, the solutions at the end of the ith generation, retrieve xlist[[c(1L, i)]]. This will be a matrix of size length(algo$min) times algo$nP. (To be consistent with other functions, xlist is itself a list. In the case of DEopt, it contains just one element.)

classify

Logical; default is FALSE. If TRUE, the result will have a class attribute TAopt 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 single-column 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 algo$storeF is TRUE, a matrix of size algo$nG times algo$nP containing the objective function values of all solutions over the generations; else NA.
xlist	if algo$storeSolutions is TRUE, a list that contains a list P of matrices and a matrix initP (the initial solution); else NA.
initial.state	the value of .Random.seed when the function was called.

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/C2017-0-01621-X

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'. https://enricoschumann.net/R/remarks.htm

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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

GAopt, PSopt

Examples

## Example 1: Trefethen's 100-digit challenge (problem 4)
## https://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/C2017-0-01621-X

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 TRUE, maximum drawdown is chosen according to percentage losses; else in units of v
summary	if TRUE, provide maximum drawdown and time when it occured; else return drawdown vector

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 v
high.position	position of high
low	the min of v
low.position	position of low

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

drawdowns

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	risk-free 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/C2017-0-01621-X

M. Gilli and Schumann, E. (2009) Implementing Binomial Trees. COMISEF Working Paper Series No. 008. https://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). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

callHestoncf

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 <- 1e-8
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 = "F-F_Research_Data_Factors_CSV.zip",
       weighting = "value", frequency = "monthly",
       price.series = FALSE, na.rm = FALSE,
       adjust.frequency = TRUE,
       return.class = "data.frame")

Arguments

dest.dir	string: a path to a directory
dataset	a character string: the CSV file name. Also supported are the keywords ‘market’ and ‘rf’.
weighting	a character string: "equal" or "value"
frequency	a character string: daily, monthly or annual. Whether it is used or ignored depends on the particular dataset.
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 TRUE, frequency is switched to ‘"daily"’ when the word ‘"daily"’ appears in the dataset's name
return.class	a string: ‘data.frame’ (the default) or ‘zoo’

Details

The function downloads data provided by Kenneth French at https://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).

For dataset ‘market’, the function downloads the three-factor dataset, and adds excess return and risk-free rate. For dataset ‘rf’, only the risk-free rate is returned.

Value

A data.frame, with contents depending on the particular dataset, or an object as specified by return.class.

If the download fails, the function evaluates to NULL, typically with warnings.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

Shiller

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, "F-F_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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 so-called population), all of which are coded as vectors containing only zeros and ones. (In GAopt, solutions are of mode logical.) The algorithm starts with randomly-chosen or user-supplied 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 randomly-changed 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) times algo$nP, or a function that creates such a matrix. If a function, it must take no arguments. If mode(mP) is not logical, then storage.mode(mP) will be tried (and a warning will be issued).

loopOF

logical. Should the OF be evaluated through a loop? Defaults to TRUE.

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 or multicore. Setting vectorised is equivalent to having algo$loopOF set to FALSE (and methodOF overrides loopOF). snow and multicore use functions clusterApply and mclapply, respectively. For snow, an object algo$cl needs to be specified (see below). For multicore, optional arguments can be passed through algo$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), a txtProgressBar is printed.

storeF

If TRUE (the default), the objective function values for every solution in every generation are stored and returned as matrix Fmat.

storeSolutions

If TRUE, the solutions (ie, binary strings) in every generation are stored and returned as a list P in list xlist (see Value section below). To check, for instance, the solutions at the end of the ith generation, retrieve xlist[[c(1L, i)]]. This will be a matrix of size algo$nB times algo$nP.

classify

Logical; default is FALSE. If TRUE, the result will have a class attribute TAopt 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 algo$storeF is TRUE, a matrix of size algo$nG times algo$nP containing the objective function values of all solutions over the generations; else NA
xlist	if algo$storeSolutions is TRUE, a list that contains a list P of matrices and a matrix initP (the initial solution); else NA.
initial.state	the value of .Random.seed when the function was called.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

DEopt, PSopt

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; ... arguments are also passed.
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. Fmat[ ,1L] contains the proposed solution over all iterations; Fmat[ ,2L] contains the accepted solutions.
xlist	a list
initial.state	the value of .Random.seed when the function was called.
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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 fun(x, ...), with x being a numeric vector or a list
levels	a list of levels for the arguments (see Examples)
...	objects passed to fun
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 lower and upper are to be expanded; see Details. Ignored when levels are explicitly specified, or when lower/upper are used and at least one has length greater than one. See Examples.
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 loop (the default), multicore or snow. See Details.
mc.control	a list containing settings that will be passed to mclapply if method is multicore. Must be a list of named elements; see the documentation of mclapply in parallel.
cl	default is NULL. If method snow is used, this must be a cluster object or an integer (the number of cores).
keepNames	logical: should the names of levels be kept?
asList	does fun expect a list? Default is FALSE.

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 fun.
minlevels	the levels that give this minimum.
values	a list. All the function values of fun.
levels	a list. All the levels for which fun was evaluated.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 (auto-expanded) and n
lower <- c(1,1); upper <- 3; n <- 4
res <- gridSearch(testFun, lower = lower, upper = upper, n = n)
res$minfun
sol$minlevels



## non-numeric 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

---

Local-Search 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

LSopt, TA.info

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; ... arguments are also passed.
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 override nS. The option is provided to makes it easier to compare and switch between functions LSopt, TAopt and SAopt.

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 when LSopt 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 solution x; it must return a changed solution.

printDetail

If TRUE (the default), information is printed. If an integer i greater then one, information is printed at very ith step.

printBar

If TRUE (the default), a txtProgressBar (from package utils) is printed). The progress bar is not shown if printDetail 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 matrix Fmat.

storeSolutions

default is FALSE. If TRUE, the solutions (ie, decision variables) in every generation are stored and returned in list xlist (see Value section below). To check, for instance, the current solution at the end of the ith generation, retrieve xlist[[c(2L, i)]].

OF.target

Numeric; when specified, the algorithm will stop when an objective-function value as low as OF.target (or lower) is achieved. This is useful when an optimal objective-function 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. Fmat[ ,1L] contains the proposed solution over all iterations; Fmat[ ,2L] contains the accepted solutions.
xlist	if algo$storeSolutions is TRUE, a list; else NA. Contains the neighbour solutions at a given iteration (xn) and the current solutions (xc). Example: Fmat[i, 2L] is the objective function value associated with xlist[[c(2L, i)]].
initial.state	the value of .Random.seed when the function was called.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 OF-value 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 OF-value 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 OF-value 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(1e-4,
              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 y (see Examples).
pad	Defaults to NULL. If not NULL, all elements of the returned moving average with position smaller than order are replaced by the value of pad. Sensible values may be NA or 0.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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")

---

Maximum-Sharpe-Ratio/Tangency Portfolio

Description

Compute maximum Sharpe-ratio 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 "qp" is supported.
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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Schumann, E. (2012) Computing the global minimum-variance portfolio. https://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 Monte-Carlo 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 exp of the final path, or return log values?
antithetic	logical: if TRUE, random numbers for only npaths/2 are drawn, and the random numbers are mirrored

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*(1-t/T)+ST*(t/T).

Value

A matrix of sample paths; each column contains the price path of an asset. Even with only a single time-step, 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

vanillaOptionEuropean

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 <- 1e-6                 ## 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 Conditional-Value-at-Risk (CVaR) Portfolios

Description

Compute minimum-CVaR 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 Value-at-Risk 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 is specified.
min.return	minimal required return. If m is not specified, the column means of R are used.
method	character. Currently, only "Rglpk" is supported.
groups	a list of group definitions
groups.wmin	a numeric vector
groups.wmax	a numeric vector
Rglpk.control	a list: settings passed to Rglpk_solve_LP

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/C2017-0-01621-X

Rockafellar, R. T. and Uryasev, S. (2000). Optimization of Conditional Value-at-Risk. Journal of Risk. 2 (3), 21–41.

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Schumann, E. (2020) Minimising Conditional Value-at-Risk (CVaR). https://enricoschumann.net/notes/minimising-conditional-var.html

See Also

minvar

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–Absolute-Deviation Portfolios

Description

Compute minimum mean–absolute-deviation 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 NULL
m	a vector of expected returns. If NULL, but min.return is not NULL, then column means are used as expected returns.
demean	logical. If TRUE, the columns of R are demeaned, corresponding to an objective function xxxx
method	string. Supported are lp and ls.
groups	group definitions
groups.wmin	list of vectors
groups.wmax	list of vectors
Rglpk.control	a list

Details

Compute the minimum mean–absolute-deviation 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) Mean-Absolute Deviation Portfolio Optimization Model and Its Applications to Tokyo Stock Market. Management Science. 37 (5), 519–531.

See Also

minvar, minCVaR

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))

---

Minimum-Variance Portfolios

Description

Compute minimum-variance 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 "qp" is supported.
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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Schumann, E. (2012) Computing the global minimum-variance portfolio. https://enricoschumann.net/R/minvar.htm

See Also

TAopt

Examples

## variance-covariance matrix from daily returns, 1 Jan 2014 -- 31 Dec 2013, of
## cleaned data set at https://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	risk-free rate
lambda	risk–reward trade-off
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 rf is specified, an additional row is added, providing the cash weight.

The i-th portfolio on the frontier corresponds to the i-th elements of return and volatility, and the i-th 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

See Also

minvar for computing the minimum-variance 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")

---

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: $\beta_1$, $\beta_2$, $\beta_3$, $\lambda$. For NSS: a vector: $\beta_1$, $\beta_2$, $\beta_3$, $\beta_4$, $\lambda_1$, $\lambda_2$.
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 Nelson-Siegel-Svensson model, COMISEF Working Paper Series No. 031. https://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/C2017-0-01621-X

Gilli, M. and Schumann, E. (2010) A Note on ‘Good’ Starting Values in Numerical Optimisation, COMISEF Working Paper Series No. 044. https://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). https://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

NSf, NSSf

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")

---

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 $\lambda$ parameter of the NS model (a scalar)
lambda1	the $\lambda_1$ parameter of the NSS model (a scalar)
lambda2	the $\lambda_2$ parameter of the NSS model (a scalar)
tm	a numeric vector with times-to-payment/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 Nelson-Siegel-Svensson model, COMISEF Working Paper Series No. 031. https://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/C2017-0-01621-X

Gilli, M. and Schumann, E. (2010) A Note on ‘Good’ Starting Values in Numerical Optimisation, COMISEF Working Paper Series No. 044. https://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). https://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

NS, NSS

Examples

## Nelson-Siegel
cor(NSf(lambda = 6, tm = 1:10)[-1L, -1L])

## Nelson-Siegel-Svensson
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 2012-02-10.

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/ex-en/ ; Euribor rates and the parameters of the Nelson-Siegel-Svensson 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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://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 column-wise.

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/C2017-0-01621-X

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Examples

pm(x <- rnorm(100), 2)
var(x)/2

pm(x, 2, normalise = TRUE)
sqrt(var(x)/2)

---