Description Details Author(s) References Examples
Functions, data and other R code from the book ‘Numerical Methods and Optimization in Finance’. Many functions are still experimental; interfaces may change in the future; there is little error handling. Comments/corrections/remarks/suggestions are very welcome (please contact the maintainer directly).
The package contains implementations of several optimisation
heuristics: Differential Evolution (DEopt
), Genetic
Algorithms (GAopt
), (Stochastic) Local Search
(LSopt
), Particle Swarm (PSopt
) and
Threshold Accepting (TAopt
). The term heuristic is
understood here as a general purpose optimisation method.
There are also various functions for option pricing; for instance,
vanillaOptionEuropean
or
vanillaOptionAmerican
.
Dependencies: The package is completely written in R. A number of
packages are suggested, but they are not necessary to use the
NMOF package. Package MASS is required to run the complete
example for PSopt
and also in one of the vignettes
(PSlms
); packages snow and parallel are optional
for functions bracketing
, GAopt
,
gridSearch
and restartOpt
, and they may
become an option for other functions; quadprog is needed for a
vignette (TAportfolio
) and some tests; RUnit is needed to
run the tests in subdirectory ‘unitTests
’.
Version numbering: minor version numbers (like in x.1-x
) are
incremented when a feature is added or an existing feature is
substantially revised. The patch level (x.x-1
) is incremented
with any published change.
Enrico Schumann
Maintainer: Enrico Schumann <es@enricoschumann.net>
Gilli, M., Maringer, D. and Schumann, E. (2011) Numerical Methods and Optimization in Finance. Elsevier. http://www.elsevierdirect.com/product.jsp?isbn=9780123756626
Schumann, E. (2013) The NMOF Manual. http://enricoschumann.net/NMOF.htm
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | ## Not run:
require("NMOF")
## overview
packageDescription("NMOF")
help(package = "NMOF")
## code from the book
showExample("equations.R")
showExample("exampleLS.R", chapter = 13)
## show NEWS file
news(Version >= "0.27-0", package = "NMOF")
## vignettes
vignette(package = "NMOF")
nss <- vignette("DEnss", package = "NMOF")
print(nss)
edit(nss)
## book website
browseURL("http://nmof.net")
browseURL("http://enricoschumann.net/NMOF.htm")
## more examples
file.show(system.file("NMOFex/README", package = "NMOF"))
file.show(system.file("NMOFex/NMOFman.R", package = "NMOF"))
## unit tests
file.show(system.file("unitTests/report.txt", package = "NMOF"))
test.rep <- readLines(system.file("unitTests/report.txt", package = "NMOF"))
nt <- gsub(".*\\((\\d+) checks?\\).*", "\\1",
grep("\\(\\d+ checks?\\)", test.rep))
message("Number of unit tests: ", sum(as.numeric(nt)))
## End(Not run)
|
Package: NMOF
Type: Package
Title: Numerical Methods and Optimization in Finance
Version: 0.40-0
Date: 2016-10-20
Maintainer: Enrico Schumann <es@enricoschumann.net>
Authors@R: person("Enrico", "Schumann", role = c("aut", "cre"), email =
"es@enricoschumann.net")
Depends: R (>= 2.14)
Imports: grDevices, graphics, parallel, stats, utils
Suggests: MASS, RUnit, quadprog
Description: Functions, examples and data from the book "Numerical
Methods and Optimization in Finance" by M. 'Gilli', D.
'Maringer' and E. Schumann (2011), ISBN 978-0123756626. The
package provides implementations of several optimisation
heuristics, such as Differential Evolution, Genetic Algorithms
and Threshold Accepting. There are also functions for the
valuation of financial instruments, such as bonds and options,
and functions that help with stochastic simulations.
License: GPL-3
URL: http://nmof.net, http://enricoschumann.net/NMOF.htm
LazyLoad: yes
LazyData: yes
Classification/JEL: C61, C63
NeedsCompilation: no
Packaged: 2016-10-20 09:47:40 UTC; es
Author: Enrico Schumann [aut, cre]
Repository: CRAN
Date/Publication: 2016-10-20 22:12:46
Built: R 3.3.1; ; "Sat, 22 Oct 2016 18:24:00 +0000"; unix
-- File: /usr/lib/R/site-library/NMOF/Meta/package.rds
Information on package 'NMOF'
Description:
Package: NMOF
Type: Package
Title: Numerical Methods and Optimization in Finance
Version: 0.40-0
Date: 2016-10-20
Maintainer: Enrico Schumann <es@enricoschumann.net>
Authors@R: person("Enrico", "Schumann", role = c("aut",
"cre"), email = "es@enricoschumann.net")
Depends: R (>= 2.14)
Imports: grDevices, graphics, parallel, stats, utils
Suggests: MASS, RUnit, quadprog
Description: Functions, examples and data from the book
"Numerical Methods and Optimization in Finance"
by M. 'Gilli', D. 'Maringer' and E. Schumann
(2011), ISBN 978-0123756626. The package
provides implementations of several optimisation
heuristics, such as Differential Evolution,
Genetic Algorithms and Threshold Accepting.
There are also functions for the valuation of
financial instruments, such as bonds and options,
and functions that help with stochastic
simulations.
License: GPL-3
URL: http://nmof.net,
http://enricoschumann.net/NMOF.htm
LazyLoad: yes
LazyData: yes
Classification/JEL: C61, C63
NeedsCompilation: no
Packaged: 2016-10-20 09:47:40 UTC; es
Author: Enrico Schumann [aut, cre]
Repository: CRAN
Date/Publication: 2016-10-20 22:12:46
Built: R 3.3.1; ; "Sat, 22 Oct 2016 18:24:00 +0000";
unix
Index:
DEopt Optimisation with Differential Evolution
EuropeanCall Computing Prices of European Calls with a
Binomial Tree
GAopt Optimisation with a Genetic Algorithm
LS.info Local-Search Information
LSopt Stochastic Local Search
MA Simple Moving Average
NMOF-package Numerical Methods and Optimization in Finance
NS Zero Rates for Nelson-Siegel-Svensson Model
NSf Factor Loadings for Nelson-Siegel and
Nelson-Siegel-Svensson
PSopt Particle Swarm Optimisation
TA.info Threshold-Accepting Information
TAopt Optimisation with Threshold Accepting
bracketing Zero-Bracketing
bundData German Government Bond Data
callCF Price a Plain-Vanilla Call with the
Characteristic Function
callHestoncf Price of a European Call under the Heston Model
callMerton Price of a European Call under Merton's
Jump-Diffusion Model
colSubset Full-rank Column Subset
drawdown Drawdown
fundData Mutual Fund Returns
gridSearch Grid Search
mc Option Pricing via Monte-Carlo Simulation
optionData Option Data
pm Partial Moments
putCallParity Put-Call Parity
qTable Prepare LaTeX Table with Quartile Plots
repairMatrix Repair an Indefinite Correlation Matrix
resampleC Resample with Specified Rank Correlation
restartOpt Restart an Optimisation Algorithm
showExample Display examples
testFunctions Classical Test Functions for Unconstrained
Optimisation
vanillaBond Pricing Plain-Vanilla Bonds
vanillaOptionEuropean Pricing Plain-Vanilla Options (European and
American)
xtContractValue Contract Value of Australian Government Bond
Future
xwGauss Integration of Gauss-type
Further information is available in the following vignettes in
directory '/usr/lib/R/site-library/NMOF/doc':
An_overview: An Overview of the NMOF Package (source, pdf)
DEnss: Fitting the Nelson--Siegel--Svensson model with Differential
Evolution (source, pdf)
LSselect: Asset selection with Local Search (source, pdf)
PSlms: Robust Regression with Particle Swarm Optimisation and
Differential Evolution (source, pdf)
TAportfolio: Portfolio Optimisation with Threshold Accepting (source,
pdf)
qTableEx: Examples for the qTable function (source, pdf)
repair: Repairing solutions (source, pdf)
vectorise: Vectorised objective functions (source, pdf)
equations.R
# equation.R -- version 2011-01-03
m <- 10000; n <- 10; trials <- 200
X <- array(rnorm(m * n), dim = c(m, n))
y <- rnorm(m)
# QR
system.time({
for (r in seq(trials)) sol1 <- qr.solve(X, y)
})
# form (X'X) and (X'y)
system.time({
for (r in seq(trials))
sol2 <- solve(crossprod(X), crossprod(X, y))
})
# cholesky
system.time({
for (r in seq(trials)) {
C <- chol(crossprod(X))
rhs <- crossprod(X, y)
sol3 <- backsolve(C, forwardsolve(t(C), rhs))
}
})
# check
stopifnot(all.equal(as.numeric(sol1), as.numeric(sol2)))
stopifnot(all.equal(as.numeric(sol2), as.numeric(sol3)))
Chapter File
1 Introduction equations.R
exampleLS.R
# exampleLS.R -- version 2010-12-21
require("NMOF")
# create random data
na <- 500L
C <- array(0.6, dim = c(na,na)); diag(C) <- 1
minVol <- 0.20; maxVol <- 0.40
Vols <- (maxVol - minVol) * runif(na) + minVol
Sigma <- outer(Vols,Vols) * C
# objective function
OF <- function(x, data) {
xx <- as.logical(x)
w <- x/sum(x)
res <- crossprod(w[xx],data$Sigma[xx,xx])
res <- tcrossprod(w[xx],res)
res
}
# neighborhood function
neighbour <- function(xc, data) {
xn <- xc
p <- sample.int(data$na, data$nn, replace = FALSE)
xn[p] <- abs(xn[p]-1L)
# reject infeasible solution
if( (sum(xn)>data$Ksup) || (sum(xn)<data$Kinf) ) {
return(xc)} else return(xn)
}
# data
data <- list(Sigma = Sigma,
Kinf = 30L,
Ksup = 60L,
na = na,
nn = 1L)
# random solution
card0 <- sample(data$Kinf:data$Ksup, 1L, replace = FALSE)
assets <- sample.int(na, card0, replace = FALSE)
x0 <- numeric(na)
x0[assets] <- 1L
# settings
algo <- list(x0 = x0,
neighbour = neighbour,
nS = 5000L)
system.time(sol1 <- LSopt(OF, algo, data))
# result
sqrt(sol1$OFvalue)
# plot best solution over time
par(ylog = TRUE); plot(sqrt(sol1$Fmat[ ,2L]), type = "l")
# run more trials
trials <- 100L
allRes <- restartOpt(LSopt, n = trials, OF, algo = algo, data = data)
allResOF <- numeric(trials)
for (i in 1L:trials) allResOF[i] <- sqrt(allRes[[i]]$OFvalue)
# a simpler objective function
OF2 <- function(x, data) {
xx <- as.logical(x); w <- 1/sum(x)
res <- sum(w * w * data$Sigma[xx,xx])
res
}
Chapter File
1 PortfolioOptimization exampleLS.R
Changes in version 0.40-0 (2016-10-20):
o new function 'xtTickValue', which computes the tick
value for of Australian government-bond futures
Changes in version 0.39-0 (2016-09-30):
o new function 'xtContractValue', for computing
contract values of Australian government-bond
futures
Changes in version 0.38-0 (2016-09-03):
o fixed: LSopt did not switch off 'printBar' in cases
when 'printDetail > 1'
o updates to documentation of LSopt/TAopt
Changes in version 0.37-6 (2016-08-27):
o fixed: 'gbb' failed when only a single path was to
be computed
Changes in version 0.37-5 (2016-08-23):
o gbm: add fanplot example
Changes in version 0.37-4 (2016-07-12):
o fixed: 'DEopt' did not work for univariate
optimisation models. (Thanks to Torsten von
Bartenwerffer for reporting.)
Changes in version 0.37-3 (2016-07-06):
o TAopt: new option 'OF.target', which lets the
algorithm stop when a desired objective-function
value is reached
Changes in version 0.37-2 (2016-06-08):
o TA.info: the best solution ('xbest') is now also
provided
Changes in version 0.37-1 (2016-05-24):
o new test function 'tfEggholder'
Changes in version 0.37-0 (2016-04-11):
o TAopt: experimental/rudimentary plot and print
methods. The methods are invoked if 'algo$classify'
is TRUE (default is FALSE in order to keep old
behaviour)
o there is now a public Git repository at
https://github.com/enricoschumann/NMOF
Changes in version 0.36-5 (2016-04-07):
o repairMatrix: new argument 'eps', which specifies
the replacement value for negative eigenvalues
Changes in version 0.36-4 (2015-11-18):
o bracketing: the resulting matrix's columns are
named 'lower' and 'upper'
Changes in version 0.36-3 (2015-11-09):
o restartOpt: when 'cl' is specified and 'method' is
identical to 'loop', the function will now issue a
warning before changing method to 'snow'
Changes in version 0.36-2 (2015-10-16):
o functions from packages other than 'base' are now
explicity imported: 'grDevices', 'graphics',
'stats', 'utils'
Changes in version 0.36-1 (2015-07-22):
o TAopt: more informative error message when
objective function evaluates to NA during
computation of thresholds
Changes in version 0.36-0 (2015-06-23):
o experimental: new function 'LS.info' for accessing
information during optimisation run (eg, current
iteration)
Changes in version 0.35-1 (2015-06-04):
o small correction in parameter description in
'callMerton' (suggested by Laurs Leth)
Changes in version 0.35-0 (2015-05-12):
o new function 'pm' for computing partial moments
Changes in version 0.34-2 (2015-05-05):
o fixed: 'DEopt' did not properly handle the
'minmaxConstr' option, which could result in
solutions that violated the box constraints.
(Thanks to David Ashear for reporting.)
Changes in version 0.34-1 (2015-02-03):
o no user-visible changes (see ChangeLog for non-visible changes)
Changes in version 0.34-0 (2014-10-23):
o Functions such as 'mclapply' from package 'parallel' are imported;
'parallel' is no longer attached when these functions are used
(ie, the package is not added to the user's search path). [Please
note that the CRAN version of NMOF has required R (>= 2.14) for
some time now, and since that R-version package 'parallel' has
been available anyway.]
o report 'Distributed computations with the NMOF package' was
updated
Changes in version 0.33-0 (2014-08-22):
o updates in documentation and vignettes
Changes in version 0.32-1 (2014-08-14):
o argument 'offset' in 'ytm' can be a vector (that worked before as
well, but is now documented)
o several small updates in documentation
Changes in version 0.32-0 (2014-03-26):
o new function 'colSubset'
o experimental: new function 'TA.info' for accessing information
during optimisation run (eg, current iteration)
o rewrote example in 'Asset selection with Local Search'
Changes in version 0.31-2 (2014-02-21):
o various small updates in documentation and vignettes (eg,
links to mailing lists)
Changes in version 0.31-1 (2013-12-31):
o expanded vignette 'Overview of the NMOF package'
Changes in version 0.31-0 (2013-12-27):
o new function 'callMerton' for European calls under Merton's
jump--diffusion model
o new vignette 'Overview of the NMOF package'
o smaller changes/updates in documentation
Changes in version 0.30-0 (2013-11-18):
o new function 'duration'
Changes in version 0.29-2 (2013-11-02):
o function 'ytm' has an argument 'offset' now; see ?ytm for examples
Changes in version 0.29-1 (2013-10-31):
o function 'ytm' is slightly faster now
o small additions to documentation
Changes in version 0.29-0 (2013-10-13):
o added function 'drawdown'
o Code examples in subdirectory 'NMOFex' (manual)
updated.
o small updates in documentation
Changes in version 0.28-2 (2013-09-04):
o vignettes are now all placed in vignettes, again [see NEWS entries
v0.27-5 (2013-04-26) and v0.26-0 (2012-09-19)]; small changes in
vignettes
Changes in version 0.28-1 (2013-07-22):
o 'snow' was removed from 'Suggested'. Whenever distribution of
computations is done via method 'snow', package 'parallel' is now
used. Note that this is currently experimental and little-tested;
feedback is thus much appreciated.
o in 'callHestoncf' and 'callCF', the upper integration limit is now
'Inf' (suggested by Wolfgang H\"ormann). This may make it easier
to compare different numerical implementations (but accuracy with
the old cutoff value was sufficient: numeric differences in prices
are of the order of a few hundredths of one cent).
Changes in version 0.28-0 (2013-06-08):
o 'multicore' was removed from the 'Suggested' packages. (Package
'parallel' had already been used -- if possible -- since version
0.25-6; see the NEWS entry.)
o the package now 'Depends' on R version >= 2.14 (because of
'parallel')
o new functions 'gbm' and 'gbb' added (geometric Brownian
motion/bridge)
Changes in version 0.27-5 (2013-04-26):
o since vignettes are built offline, the number of restarts was
increased again; cf. NEWS entry for v0.26-0 (2012-09-19)
o fixed a small error in help page for 'vanillaBond'
o internal changes in option pricing functions; more tests added for
dividend handling
Changes in version 0.27-4 (2013-04-22):
o function 'putCallParity' added
o for 'vanillaOptionEuropean', a parameter 'vol' can be passed via
'...' (which is translated into 'v' by squaring and ignored if 'v'
is already specified)
o a few internal changes and doc changes
Changes in version 0.27-3 (2013-01-04):
o Experimental: 'vanillaOptionEuropean' has arguments 'model'
(default is NULL which corresponds to Black-Scholes-Merton)
and '...'. This allows to pass other models; for
example, 'heston'. With 'greeks' set to TRUE, the function
also returns (some) derivatives. More will be added over
time, so the results should be extracted by name, not by
position (eg, x[["delta"]], but not x[[2L]]).
(Adding derivatives for non-BSM models was suggested by
Nitesh Kumar.)
Changes in version 0.27-2 (2013-01-03):
o Small corrections in docs for 'ytm' and
'vanillaOptionEuropean'.
Changes in version 0.27-1 (2013-01-02):
o Added functions 'ytm' and 'vanillaBond': compute
yield-to-maturity and price of plain-vanilla bond.
Changes in version 0.27-0 (2012-11-29):
o Function 'qTable' has new arguments 'funs', 'tabular.format'
and 'skip'; see vignette 'qTableEx' and ?qTable.
The function now (correctly) adds more space at the end of
the tabular (unless 'skip' is FALSE, which corresponds to the
previous behaviour).
Vignettes in package 'NMOF':
An_overview An Overview of the NMOF Package (source, pdf)
LSselect Asset selection with Local Search (source, pdf)
qTableEx Examples for the qTable function (source, pdf)
DEnss Fitting the Nelson--Siegel--Svensson model with
Differential Evolution (source, pdf)
TAportfolio Portfolio Optimisation with Threshold Accepting
(source, pdf)
repair Repairing solutions (source, pdf)
PSlms Robust Regression with Particle Swarm
Optimisation and Differential Evolution
(source, pdf)
vectorise Vectorised objective functions (source, pdf)
sh: 0: Can't open /dev/null
sh: 1: vi: Permission denied
sh: 0: Can't open /dev/null
sh: 0: Can't open /dev/null
The file NMOFman.R contains the code from the
NMOF manual.
The file NMOFex.R contains the code from the report
'Examples and Extensions for the NMOF package' (but
this report was superseded by the NMOF manual).
The file NMOFdist.R contains the code from the report
'Distributed computations with the NMOF package'.
The pdf-versions of the reports are available from
http://enricoschumann.net/files/NMOFman.pdf
http://enricoschumann.net/files/NMOFdist_Windows.pdf
http://enricoschumann.net/files/NMOFdist_Linux.pdf
### R code from vignette source '/home/es/Documents/Books/NMOFman/NMOFman.Rnw'
### Encoding: UTF-8
###################################################
### code chunk number 1: NMOFman.Rnw:155-157
###################################################
version <- as.Date("2016-04-20")
options(continue = " ", digits = 3, width = 60, useFancyQuotes = FALSE)
###################################################
### code chunk number 2: NMOFman.Rnw:247-249 (eval = FALSE)
###################################################
## code <- system.file("NMOFex/NMOFman.R", package = "NMOF")
## file.show(code, title = "NMOF manual")
###################################################
### code chunk number 3: NMOFman.Rnw:272-274
###################################################
require("NMOF")
set.seed(123321)
###################################################
### code chunk number 4: NMOFman.Rnw:303-304
###################################################
nrows <- 100L
###################################################
### code chunk number 5: NMOFman.Rnw:308-313
###################################################
rho <- 0.6 ## correlation between a and b
C <- array(rho, dim = c(2L, 2L))
diag(C) <- 1
ab <- array(rnorm(nrows * 2), dim = c(nrows, 2L)) %*% chol(C)
colnames(ab) <- c("a", "b")
###################################################
### code chunk number 6: NMOFman.Rnw:317-327
###################################################
par(mfrow = c(1, 3), bty = "n", mar = c(4, 4, 0, 0),
tck = 0.01, las = 1, ps = 9, mgp = c(1.8, 0.5, 0),
cex.axis = 1.2, cex.lab = 1.2, lab = c(3, 3, 7))
prx <- pretty(ab)
lims <- c(min(prx), max(prx))
plot(ab, pch = 19, cex = 0.5, asp = 1, xlim = lims, ylim = lims,
xaxt = "n", yaxt = "n", xlab = "a", ylab = "")
axis(1, col = grey(0.5))
axis(2, col = grey(0.5))
mtext("b", side = 2, line = 1.8, cex = 0.8)
###################################################
### code chunk number 7: NMOFman.Rnw:342-383
###################################################
plot_subsets <- function(ab, subset1, subset2, linesAtZero = FALSE) {
if (missing(subset2))
subset2 <- !subset1
if (cor(ab[subset1, ])[1L,2L] < cor(ab[subset2, ])[1L,2L]) {
tmp <- subset1
subset1 <- subset2
subset2 <- tmp
}
par(mfrow = c(1L, 3L), bty = "n", mar = c(4, 4, 0, 0),
tck = 0.01, las = 1, ps = 9, mgp = c(1.8,0.5,0),
cex.axis = 1.2, cex.lab = 1.2, lab = c(3,3,7))
prx <- pretty(ab)
lims <- c(min(prx), max(prx))
plot(ab[subset1, ], xlim = lims, ylim = lims, xlab = "a", ylab = "",
col = grey(0.2), pch = 19, cex = 0.5, xaxt = "n", yaxt = "n", asp = 1)
axis(1, col = grey(0.6))
lines(ab[subset2, ], col = grey(0.75), pch = 19, type = "p", cex = 0.5, asp = 1)
axis(2, col = grey(0.6))
mtext("b", side = 2, line = 1.8, cex = 0.8)
if (linesAtZero)
abline(v = 0, h = 0, col=grey(0.6))
par(mar = c(4,2,0,2))
plot(ab[subset1, ], xlim=lims, ylim = lims, xlab = "a", ylab = "",
col = grey(0.2), pch = 19, cex = 0.5, xaxt = "n", yaxt = "n", asp = 1)
axis(1, col = grey(0.6))
if (linesAtZero)
abline(v=0,h=0, col=grey(0.5))
par(mar = c(4,0,0,4))
plot(ab[subset2, ], xlim=lims, ylim = lims, xlab = "a", ylab = "",
col = grey(0.75), pch = 19, cex = 0.5, xaxt = "n", yaxt = "n", asp = 1)
axis(1, col = grey(0.6))
if (linesAtZero)
abline(v=0,h=0, col=grey(0.5))
invisible(NULL)
}
plot_subsets(ab, 1:50, 51:100)
###################################################
### code chunk number 8: NMOFman.Rnw:390-392
###################################################
cor(ab[ 1: 50, ])[1,2]
cor(ab[51:100, ])[1,2]
###################################################
### code chunk number 9: NMOFman.Rnw:415-416
###################################################
x0 <- rep(c(TRUE, FALSE), each = 50)
###################################################
### code chunk number 10: NMOFman.Rnw:420-421
###################################################
x0 <- runif(nrows) > 0.5
###################################################
### code chunk number 11: NMOFman.Rnw:424-425 (eval = FALSE)
###################################################
## ab[x, ]
###################################################
### code chunk number 12: NMOFman.Rnw:428-429 (eval = FALSE)
###################################################
## ab[!x, ]
###################################################
### code chunk number 13: NMOFman.Rnw:448-450
###################################################
OF <- function(x, ab)
-abs(cor(ab[x, ])[1L, 2L] - cor(ab[!x, ])[1L, 2L])
###################################################
### code chunk number 14: NMOFman.Rnw:454-457
###################################################
x0 <- rep(c(TRUE, FALSE), each = 50L)
OF( x0, ab)
OF(!x0, ab) ## should give the same result
###################################################
### code chunk number 15: NMOFman.Rnw:476-488
###################################################
trials <- 1e5
OFvalues <- numeric(trials)
solutions <- vector("list", trials)
for (i in seq_len(trials)) {
c1 <- sample(21:80, 1L) ## cardinality of subset 1
x0 <- logical(nrows)
x0[sample.int(nrows, c1)] <- TRUE
OFvalues[i] <- OF(x0, ab) ## store results
solutions[[i]] <- x0
}
###################################################
### code chunk number 16: NMOFman.Rnw:492-493
###################################################
summary(OFvalues)
###################################################
### code chunk number 17: NMOFman.Rnw:497-499
###################################################
xbest <- which.min(OFvalues)
OFvalues[xbest]
###################################################
### code chunk number 18: NMOFman.Rnw:506-507
###################################################
xRandom <- solutions[[xbest]]
###################################################
### code chunk number 19: NMOFman.Rnw:524-526
###################################################
subset1 <- ab[ ,1L] * ab[ ,2L] > 0
subset2 <- ab[ ,1L] * ab[ ,2L] <= 0
###################################################
### code chunk number 20: NMOFman.Rnw:531-532
###################################################
plot_subsets(ab, subset1, linesAtZero=TRUE)
###################################################
### code chunk number 21: NMOFman.Rnw:545-546
###################################################
OF(subset1, ab)
###################################################
### code chunk number 22: NMOFman.Rnw:551-553
###################################################
sum(subset1)
sum(subset2)
###################################################
### code chunk number 23: NMOFman.Rnw:560-567
###################################################
cr <- order(ab[ ,1L] * ab[ ,2L])
OFvalues <- numeric(nrows)
for (i in 20:80) {
x0 <- logical(nrows)
x0[cr[seq_len(i)]] <- TRUE
OFvalues[i] <- OF(x0, ab)
}
###################################################
### code chunk number 24: NMOFman.Rnw:570-575
###################################################
cutoff <- which.min(OFvalues)
subset1 <- logical(nrows)
subset1[cr[seq_len(nrows) <= cutoff]] <- TRUE
subset2 <- !subset1
OF(subset1, ab)
###################################################
### code chunk number 25: NMOFman.Rnw:580-581
###################################################
xConstr <- subset1
###################################################
### code chunk number 26: NMOFman.Rnw:618-653
###################################################
greedy <- function(fun, x0, ab, n, nmin, maxit = 1000L) {
done <- FALSE
xbest <- xc <- x0
xbestF <- xcF <- fun(xbest, ab)
ic <- 0
while (!done) {
if (ic > maxit)
break
else
ic <- ic + 1L
done <- TRUE
xc <- xbest
for (i in seq_len(n)) {
xn <- xc ## create a new solution
xn[i] <- !xn[i]
sxn <- sum(xn) ## check constraints
enough <- sxn >= nmin
notTooMany <- sxn <= n - nmin
if (enough && notTooMany) {
xnF <- fun(xn, ab)
if (xnF < xbestF) {
xbest <- xn
xbestF <- xnF
done <- FALSE
}
}
}
}
list(xbest = xbest, OFvalue = xbestF, ic = ic)
}
###################################################
### code chunk number 27: NMOFman.Rnw:660-667
###################################################
x0 <- rep(c(TRUE, FALSE), each = 50L)
result <- greedy(fun = OF, x0 = x0, ab = ab,
n = 100, nmin = 20, maxit = 1000L)
xGreedy <- result$xbest
OF(x0, ab)
OF(xGreedy, ab)
result$OFvalue
###################################################
### code chunk number 28: NMOFman.Rnw:670-671
###################################################
result$ic
###################################################
### code chunk number 29: NMOFman.Rnw:677-691
###################################################
trials <- 1000
OFvalues <- numeric(trials)
solutions <- vector("list", trials)
moves <- numeric(trials)
for (i in seq_len(trials)) {
x0 <- runif(nrows) > 0.5
result <- greedy(fun = OF, x0 = x0, ab = ab,
n = 100, nmin = 20, maxit = 1000L)
OFvalues[i] <- result$OFvalue
solutions[[i]] <- result$xbest
moves[i] <- result$ic
}
###################################################
### code chunk number 30: NMOFman.Rnw:695-696
###################################################
summary(OFvalues)
###################################################
### code chunk number 31: NMOFman.Rnw:703-704
###################################################
summary(moves)
###################################################
### code chunk number 32: NMOFman.Rnw:712-715
###################################################
xbest <- which.min(OFvalues)
OFvalues[xbest]
xGreedy <- solutions[[xbest]]
###################################################
### code chunk number 33: NMOFman.Rnw:739-740
###################################################
Data <- list(ab = ab, nrows = nrows, nmin = 20L)
###################################################
### code chunk number 34: NMOFman.Rnw:745-746
###################################################
x0 <- runif(nrows) > 0.5
###################################################
### code chunk number 35: NMOFman.Rnw:753-768
###################################################
neighbour <- function(xc, Data) {
xn <- xc
p <- sample.int(Data$nrows, size = 1L)
xn[p] <- !xn[p]
sxn <- sum(xn)
enough <- sxn >= Data$nmin
notTooMany <- sxn <= (Data$nrows - Data$nmin)
if (enough && notTooMany)
xn
else
xc
}
###################################################
### code chunk number 36: NMOFman.Rnw:777-778
###################################################
table(x0 == neighbour(x0, Data))
###################################################
### code chunk number 37: NMOFman.Rnw:781-783
###################################################
OF <- function(x, Data)
-abs(cor(Data$ab[x, ])[1L, 2L] - cor(Data$ab[!x, ])[1L, 2L])
###################################################
### code chunk number 38: NMOFman.Rnw:786-788
###################################################
OF(x0, Data)
OF(neighbour(x0, Data), Data)
###################################################
### code chunk number 39: NMOFman.Rnw:792-798
###################################################
algo <- list(nS = 3000L, ## number of steps to make
neighbour = neighbour, ## neighbourhood function
x0 = x0, ## initial solution
printBar = FALSE)
sol1 <- LSopt(OF, algo = algo, Data = Data)
(sol1$OFvalue)
###################################################
### code chunk number 40: NMOFman.Rnw:803-807
###################################################
subset1 <- sol1$xbest
subset2 <- !sol1$xbest
c1 <- cor(ab[subset1, ])[1L, 2L]
c2 <- cor(ab[subset2, ])[1L, 2L]
###################################################
### code chunk number 41: NMOFman.Rnw:811-812
###################################################
plot_subsets(ab, subset1)
###################################################
### code chunk number 42: ranfun
###################################################
makex0 <- function() {
x <- logical(100)
while (sum(x) > 80 || sum(x) < 20)
x <- runif(100) > 0.5
x
}
###################################################
### code chunk number 43: NMOFman.Rnw:837-848
###################################################
algo <- list(nS = 4000L, ## number of steps to make
neighbour = neighbour, ## neighbourhood function
x0 = makex0, ## initial solution
printBar = FALSE,
printDetail = FALSE)
restarts1 <- restartOpt(LSopt, 100, OF = OF, algo = algo, Data)
restarts1OFvalues <- sapply(restarts1, `[[`, "OFvalue")
algo$nS <- 8000L
restarts2 <- restartOpt(LSopt, 100, OF = OF, algo = algo, Data)
restarts2OFvalues <- sapply(restarts2, `[[`, "OFvalue")
###################################################
### code chunk number 44: NMOFman.Rnw:853-859
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0), ps = 8, tck = 0.001)
plot( ecdf(restarts1OFvalues), main = "", ylab = "", xlab = "",
cex = 0.4, pch = 19, col = grey(.2), xlim = c(-2,-1))
lines(ecdf(restarts2OFvalues),
cex = 0.4, pch = 19, col = grey(.6))
abline(v = OF(xConstr, Data))
###################################################
### code chunk number 45: NMOFman.Rnw:879-884
###################################################
algo$q <- 0.9
algo$nT <- 10
algo$nS <- 400
sol2 <- TAopt(OF, algo = algo, Data = Data)
sol2$OFvalue
###################################################
### code chunk number 46: NMOFman.Rnw:888-893
###################################################
subset1 <- sol2$xbest
subset2 <- !sol2$xbest
c1 <- cor(ab[subset1, ])[1L, 2L]
c2 <- cor(ab[subset2, ])[1L, 2L]
OF(sol2$xbest, Data)
###################################################
### code chunk number 47: NMOFman.Rnw:897-898
###################################################
plot_subsets(ab, subset1)
###################################################
### code chunk number 48: NMOFman.Rnw:903-911
###################################################
algo$printBar <- FALSE
algo$printDetail <- FALSE
restarts3 <- restartOpt(TAopt, 100, OF = OF, algo = algo, Data)
restarts3OFvalues <- sapply(restarts3, `[[`, "OFvalue")
algo$nS <- 1500
restarts4 <- restartOpt(TAopt, 100, OF = OF, algo = algo, Data)
restarts4OFvalues <- sapply(restarts4, `[[`, "OFvalue")
###################################################
### code chunk number 49: NMOFman.Rnw:914-924
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0), ps = 8, tck = 0.001)
plot( ecdf(restarts1OFvalues), main = "", ylab = "", xlab = "",
cex = 0.4, pch = 19, col = grey(.2), xlim = c(-2,-1))
lines(ecdf(restarts2OFvalues),
cex = 0.4, pch = 19, col = grey(.6))
lines(ecdf(restarts3OFvalues),
cex = 0.4, pch = 19, col = grey(.2), lty=2)
lines(ecdf(restarts4OFvalues),
cex = 0.4, pch = 19, col = grey(.6), lty=2)
abline(v = OF(xConstr, Data))
###################################################
### code chunk number 50: NMOFman.Rnw:942-957
###################################################
neighbour <- function(xc, Data) {
xn <- xc
p <- sample.int(Data$nrows, size = Data$size)
xn[p] <- !xn[p]
sxn <- sum(xn)
enough <- sxn >= Data$nmin
notTooMany <- sxn <= (Data$nrows - Data$nmin)
if (enough && notTooMany)
xn
else
xc
}
###################################################
### code chunk number 51: NMOFman.Rnw:964-984
###################################################
neighbour.maker <- function(n, nmin, size) {
force(n)
force(nmin)
function(xc) {
xn <- xc
p <- sample.int(n, size = size)
xn[p] <- !xn[p]
sxn <- sum(xn)
enough <- sxn >= nmin
notTooMany <- sxn <= (n - nmin)
if (enough && notTooMany)
xn
else
xc
}
}
###################################################
### code chunk number 52: NMOFman.Rnw:988-990
###################################################
N <- neighbour.maker(n = 10, nmin = 2, size = 3)
N
###################################################
### code chunk number 53: NMOFman.Rnw:992-1001
###################################################
compareLogicals <- function(x, y, ...) {
argsL <- list(...)
if (!("sep" %in% names(argsL)))
argsL$sep <- ""
do.call("cat",
c(list("\n",as.integer(x), "\n", as.integer(y), "\n",
ifelse(x == y, " ", "^"), "\n"), argsL))
message("The vectors differ in ", sum(x != y), " place(s).")
}
###################################################
### code chunk number 54: NMOFman.Rnw:1005-1008
###################################################
x0 <- rep(c(TRUE, FALSE), each = 5L)
Data <- list(nrows = 10, nmin = 2, size = 1)
compareLogicals(x0, neighbour(x0, Data))
###################################################
### code chunk number 55: NMOFman.Rnw:1012-1016
###################################################
N <- neighbour.maker(n = 10, nmin = 2, size = 1)
compareLogicals(x0, N(x0))
compareLogicals(x0, N(x0))
compareLogicals(x0, N(x0))
###################################################
### code chunk number 56: NMOFman.Rnw:1020-1024
###################################################
N <- neighbour.maker(n = 10, nmin = 2, size = 3)
compareLogicals(x0, N(x0))
compareLogicals(x0, N(x0))
compareLogicals(x0, N(x0))
###################################################
### code chunk number 57: NMOFman.Rnw:1038-1041
###################################################
require("NMOF")
require("rbenchmark")
set.seed(46457)
###################################################
### code chunk number 58: NMOFman.Rnw:1076-1085
###################################################
randomData <- function(p, n, rscale = 0.5) {
X <- array(rnorm(n * p), dim = c(n, p))
k <- sample.int(p, 1L) ## the number of regressors
K <- sample.int(p, k) ## the set of regressors
betatrue <- numeric(p)
betatrue[K] <- rnorm(k) ## the true coefficients
y <- X %*% betatrue + rnorm(n)*rscale
list(X = X, y = y, betatrue = betatrue, K = K, n = n, p = p)
}
###################################################
### code chunk number 59: NMOFman.Rnw:1093-1096
###################################################
n <- 60L
p <- 5L
rD <- randomData(p, n)
###################################################
### code chunk number 60: NMOFman.Rnw:1108-1109
###################################################
b0 <- rnorm(p)
###################################################
### code chunk number 61: NMOFman.Rnw:1120-1124
###################################################
Data <- list(X = rD$X,
y = rD$y,
p = rD$p,
n = rD$n)
###################################################
### code chunk number 62: NMOFman.Rnw:1130-1134
###################################################
OFls <- function(b, Data) {
tmp <- Data$y - Data$X %*% b
sum(tmp * tmp)
}
###################################################
### code chunk number 63: NMOFman.Rnw:1140-1145
###################################################
tmp <- rnorm(1e4)
benchmark(ignore1 <- tmp * tmp * tmp,
ignore2 <- tmp^3,
columns = c("test", "elapsed", "relative"),
replications = 1e3, order = "relative")
###################################################
### code chunk number 64: NMOFman.Rnw:1150-1151
###################################################
all.equal(ignore1, ignore2)
###################################################
### code chunk number 65: NMOFman.Rnw:1156-1162
###################################################
algo <- list(nG = 200, ## number of generations
nP = 50, ## population size
min = rep(-20, p),
max = rep( 20, p),
printBar = FALSE)
resDE <- DEopt(OFls, algo = algo, Data = Data)
###################################################
### code chunk number 66: NMOFman.Rnw:1168-1170
###################################################
data.frame(QR = qr.solve(Data$X, Data$y),
DE = resDE$xbest)
###################################################
### code chunk number 67: NMOFman.Rnw:1181-1184
###################################################
algo <- list(nP = 50,
min = rep(-20, p), max = rep( 20, p),
printBar = FALSE, printDetail = FALSE)
###################################################
### code chunk number 68: NMOFman.Rnw:1190-1198
###################################################
algo$nG <- 25
results1 <- restartOpt(DEopt, n = 50, OF= OFls, algo = algo, Data = Data)
algo$nG <- 50
results2 <- restartOpt(DEopt, n = 50, OF= OFls, algo = algo, Data = Data)
algo$nG <- 100
results3 <- restartOpt(DEopt, n = 50, OF= OFls, algo = algo, Data = Data)
algo$nG <- 200
results4 <- restartOpt(DEopt, n = 50, OF= OFls, algo = algo, Data = Data)
###################################################
### code chunk number 69: NMOFman.Rnw:1204-1209
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0), ps = 8, tck = 0.001)
plot( ecdf(sapply(results1, `[[`, "OFvalue")), main = "", ylab = "", xlab = "",
cex = 0.4, pch = 19, col = grey(.6), xlim = c(0,150), ylim = c(0,1))
lines(ecdf(sapply(results2, `[[`, "OFvalue")),
cex = 0.4, pch = 19, col = grey(.4))
###################################################
### code chunk number 70: NMOFman.Rnw:1214-1219
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0), ps = 8, tck = 0.001)
plot( ecdf(sapply(results2, `[[`, "OFvalue")), main = "", ylab = "", xlab = "",
cex = 0.4, pch = 19, col = grey(.6), ylim = c(0,1))
lines(ecdf(sapply(results3, `[[`, "OFvalue")),
cex = 0.4, pch = 19, col = grey(.4))
###################################################
### code chunk number 71: NMOFman.Rnw:1224-1229
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0), ps = 8, tck = 0.001)
plot( ecdf(sapply(results3, `[[`, "OFvalue")), main = "", ylab = "", xlab = "",
cex = 0.4, pch = 19, col = grey(.4), ylim = c(0,1))
lines(ecdf(sapply(results4, `[[`, "OFvalue")),
cex = 0.4, pch = 19, col = grey(.2))
###################################################
### code chunk number 72: NMOFman.Rnw:1241-1246
###################################################
OFlts <- function(b, Data) {
tmp <- Data$y - Data$X %*% b
tmp <- sort(tmp * tmp, partial = Data$h)
sum(tmp[seq_len(Data$h)])
}
###################################################
### code chunk number 73: NMOFman.Rnw:1254-1258
###################################################
require("robustbase")
alpha <- 0.9
Data <- list(X = rD$X, y = rD$y, p = rD$p, n = rD$n,
h = h.alpha.n(alpha, n = n, p = p))
###################################################
### code chunk number 74: NMOFman.Rnw:1265-1270
###################################################
resDE <- DEopt(OFlts, algo = algo, Data = Data)
resLTS <- ltsReg(rD$y ~ -1 + rD$X, alpha = alpha,
use.correction = FALSE)
data.frame(fastLTS = resLTS$raw.coefficients,
DE = resDE$xbest)
###################################################
### code chunk number 75: NMOFman.Rnw:1274-1280
###################################################
cLTS <- resLTS$raw.coefficients
cat("LTS")
sum(sort((Data$X %*%cLTS - Data$y)^2)[1:Data$h])
cDE <- resDE$xbest
cat("DEopt")
sum(sort((Data$X %*%cDE - Data$y)^2)[1:Data$h])
###################################################
### code chunk number 76: NMOFman.Rnw:1288-1290
###################################################
any(b0 < 0)
sum(b0)
###################################################
### code chunk number 77: NMOFman.Rnw:1299-1303
###################################################
repair <- function(b, Data) {
b <- abs(b)
b/sum(b)
}
###################################################
### code chunk number 78: NMOFman.Rnw:1310-1313
###################################################
b1 <- repair(b0, Data)
all(b1 >= 0) ## should be TRUE
sum(b1) ## should be 1
###################################################
### code chunk number 79: NMOFman.Rnw:1323-1326
###################################################
b0
b0 - abs(b0)
sum(b0)
###################################################
### code chunk number 80: NMOFman.Rnw:1332-1338
###################################################
b <- rnorm(1000L)
benchmark(ignore1 <- sum(b - abs(b))/2,
ignore2 <- sum(b[b < 0]),
columns = c("test", "elapsed", "relative"),
replications = 1e4, order = "relative")
all.equal(ignore1, ignore2)
###################################################
### code chunk number 81: NMOFman.Rnw:1342-1343
###################################################
abs(sum(b0) - 1)
###################################################
### code chunk number 82: NMOFman.Rnw:1348-1354
###################################################
Data$pw1 <- 100
Data$pw2 <- 100
penalty <- function(b, Data)
Data$pw1 * -sum(b - abs(b)) + Data$pw2 * abs(sum(b) - 1)
penalty(b0, Data)
penalty(b1, Data)
###################################################
### code chunk number 83: NMOFman.Rnw:1359-1363
###################################################
algo <- list(nG = 300, nP = 50,
min = rep(-20, p), max = rep( 20, p),
printBar = FALSE)
resDE <- DEopt(OFls, algo = algo, Data = Data)
###################################################
### code chunk number 84: NMOFman.Rnw:1368-1372
###################################################
round(resDE$xbest, 5)
resDE$OFvalue
sum(resDE$xbest)
all(resDE$xbest >= 0)
###################################################
### code chunk number 85: NMOFman.Rnw:1376-1382
###################################################
algo$repair <- repair
resDE <- DEopt(OFls, algo = algo, Data = Data)
round(resDE$xbest,5)
resDE$OFvalue
sum(resDE$xbest)
all(resDE$xbest >= 0)
###################################################
### code chunk number 86: NMOFman.Rnw:1385-1392
###################################################
algo$repair <- NULL
algo$pen <- penalty
resDE <- DEopt(OFls, algo = algo, Data = Data)
round(resDE$xbest,20)
resDE$OFvalue
sum(resDE$xbest)
all(resDE$xbest >= 0)
###################################################
### code chunk number 87: NMOFman.Rnw:1407-1408
###################################################
OFlts
###################################################
### code chunk number 88: NMOFman.Rnw:1411-1414
###################################################
b0 <- rnorm(p)
b1 <- rnorm(p)
P <- cbind(b0 = b0, b1 = b1)
###################################################
### code chunk number 89: NMOFman.Rnw:1418-1421
###################################################
head(Data$y - Data$X %*% b0)
head(Data$y - Data$X %*% b1)
head(drop(Data$y) - Data$X %*% P)
###################################################
### code chunk number 90: NMOFman.Rnw:1430-1433
###################################################
head(Data$y - Data$X %*% b0)^2
head(Data$y - Data$X %*% b1)^2
head((drop(Data$y) - Data$X %*% P)*(drop(Data$y) - Data$X %*% P))
###################################################
### code chunk number 91: NMOFman.Rnw:1436-1442
###################################################
OFlts2 <- function(b, Data) {
tmp <- drop(Data$y) - Data$X %*% b
tmp <- tmp * tmp
tmp <- apply(tmp, 2L, sort, partial = Data$h)
.colSums(tmp[seq_len(Data$h), ,drop = FALSE], Data$h, ncol(b))
}
###################################################
### code chunk number 92: NMOFman.Rnw:1446-1449
###################################################
nP <- 100
P <- array(rnorm(p * nP), dim = c(p, nP))
sol0 <- OFlts2(P, Data)
###################################################
### code chunk number 93: NMOFman.Rnw:1452-1459
###################################################
sol1 <- numeric(nP)
benchmark(for (i in seq_len(nP))
sol1[i] <- OFlts(P[ , i, drop = FALSE], Data),
sol2 <- OFlts2(P, Data),
columns = c("test", "elapsed", "relative"),
replications = 100, order = "relative")
all.equal(ignore1, ignore2)
###################################################
### code chunk number 94: NMOFman.Rnw:1858-1865
###################################################
OF <- tfTrefethen
n <- 100L
surf <- matrix(NA, n, n)
x1 <- seq(from = -10, to = 10, length.out = n)
for (i in seq_len(n))
for (j in seq_len(n))
surf[i, j] <- tfTrefethen(c(x1[i], x1[j]))
###################################################
### code chunk number 95: NMOFman.Rnw:1873-1879
###################################################
par(bty = "n", las = 1, mar = c(3,4,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
contour(x1, x1, surf, nlevels=5, col = grey(0.6))
## the actual minimum
abline(v = -0.02440308, h = 0.21061243, col = grey(0.6))
###################################################
### code chunk number 96: NMOFman.Rnw:1885-1896
###################################################
algo <- list(nP = 50L,
nG = 300L,
F = 0.6,
CR = 0.9,
min = c(-10,-10),
max = c(10,10),
printDetail = FALSE,
printBar = FALSE,
storeF = TRUE,
storeSolutions = TRUE)
sol <- DEopt(OF = OF, algo = algo)
###################################################
### code chunk number 97: NMOFman.Rnw:1901-1907
###################################################
names(sol)
sd(sol$popF)
ts.plot(sol$Fmat, xlab = "generations", ylab = "OF")
length(sol$xlist)
xlist <- sol$xlist[[1L]]
###################################################
### code chunk number 98: NMOFman.Rnw:1916-1931
###################################################
## show solution 1 (column 1) in population over time
xlist[[ 1L]][ ,1L] ## at the end of generation 1
## ...
xlist[[ 10L]][ ,1L] ## at the end of generation 10
## ...
xlist[[300L]][ ,1L] ## at the end of generation 300
res <- sapply(xlist, `[`, 1:2, 1) ## get row 1 and 2 from column 1
res2 <- sapply(xlist, `[`, TRUE, 1) ## simpler
all.equal(res, res2)
dim(res)
res[ ,1L]
res[ ,2L]
res[ ,300L]
###################################################
### code chunk number 99: NMOFman.Rnw:1937-1951
###################################################
## show parameter 2 (row 2) in population over time
xlist[[ 1L]][2L, ] ## at the end of generation 1
## ...
xlist[[ 10L]][2L, ] ## at the end of generation 10
## ...
xlist[[300L]][2L, ] ## at the end of generation 300
res <- sapply(xlist, `[`, 2, 1:50)
res <- sapply(xlist, `[`, 2, TRUE) ## simpler
dim(res)
res[ ,1L]
res[ ,2L]
res[ ,300L]
###################################################
### code chunk number 100: NMOFman.Rnw:1957-1966 (eval = FALSE)
###################################################
## ## transposing xlist[[i]] gives a two-column matrix -- see ?points
## ## initial solutions
## points(t(xlist[[1L]]), pch = 21, bg=grey(0.9), col = grey(.2))
##
## ## solutions at the end of generation 100
## points(t(xlist[[100L]]), pch = 21, bg=grey(0.9), col = grey(.2))
##
## ## solutions at the end of generation 100
## points(t(xlist[[300L]]), pch = 21, bg=grey(0.9), col = grey(.2))
###################################################
### code chunk number 101: NMOFman.Rnw:1970-1999
###################################################
setEPS()
postscript(file = "figures/c1.eps", width = 2, height = 2)
par(bty = "n", las = 1,mar = c(2,2,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
contour(x1, x1, surf, nlevels=5, col = grey(0.6))
abline(v = -0.02440308, h = 0.21061243, col = grey(0.6))
points(t(xlist[[1L]]), pch = 21, bg=grey(0.9), col = grey(.2))
invisible(dev.off())
setEPS()
postscript(file = "figures/c2.eps", width = 2, height = 2)
par(bty = "n", las = 1,mar = c(2,2,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
contour(x1, x1, surf, nlevels=5, col = grey(0.6))
abline(v = -0.02440308, h = 0.21061243, col = grey(0.6))
points(t(xlist[[100L]]), pch = 21, bg=grey(0.9), col = grey(.2))
invisible(dev.off())
setEPS()
postscript(file = "figures/c3.eps", width = 2, height = 2)
par(bty="n", las = 1,mar = c(2,2,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
contour(x1, x1, surf, nlevels=5, col = grey(0.6))
abline(v = -0.02440308, h = 0.21061243, col = grey(0.6))
points(t(xlist[[300L]]), pch = 21, bg=grey(0.9), col = grey(.2))
invisible(dev.off())
cat("\\includegraphics{figures/c1.eps}",
"\\includegraphics{figures/c2.eps}",
"\\includegraphics{figures/c3.eps}", sep = "")
###################################################
### code chunk number 102: NMOFman.Rnw:2013-2028
###################################################
OF <- function(par, Data) {
## compute model yields
y <- Data$model(par, Data$tm)
## all rates finite?
validRates <- !any(is.na(y))
if (validRates) {
## any rates negative? if yes, add penalty
pen1 <- sum(abs(y - abs(y))) * Data$ww
F <- max(abs(y - Data$yM)) + pen1
} else F <- 1e8
F
}
###################################################
### code chunk number 103: NMOFman.Rnw:2035-2052
###################################################
algo <- list(nP = 200L, nG = 100L,
F = 0.50, CR = 0.99,
min = c( 0,-10,-10, 0),
max = c( 1, 10, 10, 10),
storeSolutions = TRUE, printBar = FALSE)
## set up yield curve and put information in Data
tm <- 1:20 ## times to maturity
parTRUE <- c(5, 3, 2, 1) ## true parameters
yM <- NS(parTRUE, tm) ## true market yields
Data <- list(yM = yM, tm = tm, model = NS, ww = 0.1, maxb1 = 4)
## solve with DEopt
sol <- DEopt(OF = OF, algo = algo, Data = Data)
P <- sol$xlist[[1L]] ## all population matrices
p1 <- sapply(P, `[`, 1L, TRUE)
###################################################
### code chunk number 104: NMOFman.Rnw:2059-2068
###################################################
par(bty = "n", las = 1, mar = c(4,4,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
plot(jitter(rep(seq_len(algo$nG), each = algo$nP), factor = 5),
p1,
pch = 21, cex = 0.01, ylim = c(-5,10),
xlab = "", ylab = "")
mtext("generation", 1, line = 2)
mtext("parameter\nvalue", 2, line = 1)
###################################################
### code chunk number 105: NMOFman.Rnw:2076-2111
###################################################
OF2 <- function(par, Data) {
## compute model yields
y <- Data$model(par, Data$tm)
## all rates finite?
validRates <- !any(is.na(y))
if (validRates) {
## any rates negative? if yes, add penalty
pen1 <- sum(abs(y - abs(y))) * Data$ww
## is b1 greater than Data$maxb1? if yes, add penalty
pen2 <- par[1L] - Data$maxb1
pen2 <- pen2 + abs(pen2)
pen2 <- pen2
F <- max(abs(y - Data$yM)) + pen1 + pen2
} else F <- 1e8
F
}
## solve with DEopt
sol <- DEopt(OF = OF2, algo = algo, Data = Data)
P <- sol$xlist[[1L]] ### all population matrices
p1 <- sapply(P, `[`, 1, TRUE)
par(bty = "n", las = 1, mar = c(4,4,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
plot(jitter(rep(seq_len(algo$nG), each = algo$nP), factor = 5),
p1,
pch = 21, cex = 0.01, ylim = c(-5,10),
xlab = "", ylab = "")
abline(h = 4, col=grey(0.5))
mtext("generation", 1, line = 2)
mtext("parameter\nvalue", 2, line = 1)
###################################################
### code chunk number 106: NMOFman.Rnw:2202-2203
###################################################
tfRosenbrock
###################################################
### code chunk number 107: NMOFman.Rnw:2210-2214
###################################################
OF <- tfRosenbrock ## see ?testFunctions
size <- 5L ## set dimension
x <- rep.int(1, size) ## the known solution ...
OF(x) ## ... should give zero
###################################################
### code chunk number 108: NMOFman.Rnw:2222-2229
###################################################
algo <- list(printBar = FALSE,
nP = 50L,
nG = 500L,
F = 0.6,
CR = 0.9,
min = rep(-100, size),
max = rep( 100, size))
###################################################
### code chunk number 109: NMOFman.Rnw:2236-2242
###################################################
## a vectorised OF: works only with *matrix* x
OF2 <- function(x) {
n <- dim(x)[1L]
xi <- x[1L:(n - 1L), ]
colSums(100 * (x[2L:n, ] - xi * xi)^2 + (1 - xi)^2)
}
###################################################
### code chunk number 110: NMOFman.Rnw:2247-2251
###################################################
x <- matrix(rnorm(size * algo$nP), size, algo$nP)
c(OF(x[ ,1L]), OF(x[ ,2L]), OF(x[ ,3L]))
OF2(x)[1L:3L] ## should give the same result
all.equal(OF2(x)[1L:3L], c(OF(x[ ,1L]), OF(x[ ,2L]), OF(x[ ,3L])))
###################################################
### code chunk number 111: NMOFman.Rnw:2262-2268
###################################################
set.seed(1223445)
(t1 <- system.time(sol <- DEopt(OF = OF, algo = algo)))
algo$loopOF <- FALSE
set.seed(1223445)
(t2 <- system.time(sol2 <- DEopt(OF = OF2, algo = algo)))
###################################################
### code chunk number 112: NMOFman.Rnw:2273-2276
###################################################
sol$OFvalue ## both should be zero (with luck)
sol2$OFvalue
t1[[3L]]/t2[[3L]] ## speedup
###################################################
### code chunk number 113: NMOFman.Rnw:2308-2319
###################################################
na <- 100L ## number of assets
np <- 100L ## size of population
trials <- seq_len(100L) ## for speed test
## a covariance matrix
Sigma <- array(0.7, dim = c(na, na)); diag(Sigma) <- 1
## set up population
W <- array(runif(na * np), dim = c(na, np))
## budget constraint
scaleFun <- function(x) x/sum(x); W <- apply(W, 2L, scaleFun)
###################################################
### code chunk number 114: NMOFman.Rnw:2323-2343
###################################################
## variant 1
t1 <- system.time({
for (i in trials) {
res1 <- numeric(np)
for (j in seq_len(np)) {
w <- W[ ,j]
res1[j] <- w %*% Sigma %*% w
}
}
})
## variant 2
t2 <- system.time({
for (i in trials) res2 <- diag(t(W) %*% Sigma %*% W)
})
## variant 3
t3 <- system.time({
for (i in trials) res3 <- colSums(Sigma %*% W * W)
})
###################################################
### code chunk number 115: NMOFman.Rnw:2348-2350
###################################################
all.equal(res1,res2)
all.equal(res2,res3)
###################################################
### code chunk number 116: NMOFman.Rnw:2355-2358
###################################################
t1 ## speedup for variant 1
t2 ## speedup for variant 2
t3 ## speedup for variant 3
###################################################
### code chunk number 117: NMOFman.Rnw:2383-2397
###################################################
n <- 100L ## number of observation
p <- 5L ## number of regressors
np <- 100L ## population size
trials <- seq_len(1000L)
## random data
X <- array(rnorm(n * p), dim = c(n, p))
y <- rnorm(n)
## random population
Theta <- array(rnorm(p * np), dim = c(p, np))
## empty residuals matrix
R1 <- array(NA, dim = c(n, np))
###################################################
### code chunk number 118: NMOFman.Rnw:2401-2410
###################################################
system.time({
for (i in trials)
for (j in seq_len(np))
R1[ ,j] <- y - X %*% Theta[ ,j]
})
system.time({
for (i in trials)
R2 <- y - X %*% Theta
})
###################################################
### code chunk number 119: NMOFman.Rnw:2417-2418
###################################################
all.equal(R1, R2) ## ... should be TRUE
###################################################
### code chunk number 120: NMOFman.Rnw:2453-2462
###################################################
testFun <- function(x) {
Sys.sleep(0.1) ## wasting time...
cos(1/x^2)
}
system.time(sol1 <- bracketing(testFun, interval = c(0.3, 0.9),
n = 100L))
system.time(sol2 <- bracketing(testFun, interval = c(0.3, 0.9),
n = 100L, method = "snow", cl = 2))
all.equal(sol1, sol2)
###################################################
### code chunk number 121: NMOFman.Rnw:2486-2490
###################################################
testFun <- function(x) {
Sys.sleep(0.1) ## wasting time...
x[1L] + x[2L]^2
}
###################################################
### code chunk number 122: NMOFman.Rnw:2494-2498
###################################################
lower <- c(1, 3); upper <- 5; n <- 5L
system.time(sol1 <- gridSearch(fun = testFun,
lower = lower, upper = upper,
n = n, printDetail = TRUE))
###################################################
### code chunk number 123: NMOFman.Rnw:2502-2504
###################################################
seq(from = 1, to = 5, length.out= n) ## x_1
seq(from = 3, to = 5, length.out= n) ## x_2
###################################################
### code chunk number 124: NMOFman.Rnw:2508-2510
###################################################
sol1$minfun
sol1$minlevels
###################################################
### code chunk number 125: NMOFman.Rnw:2514-2521
###################################################
system.time(sol2 <- gridSearch(fun = testFun,
lower = lower,
upper = upper,
n = n, printDetail = FALSE,
method = "snow", ## use 'snow' ...
cl = 2L)) ## ... with 2 cores
all.equal(sol1, sol2)
###################################################
### code chunk number 126: NMOFman.Rnw:2561-2562
###################################################
cfBSM
###################################################
### code chunk number 127: NMOFman.Rnw:2566-2585
###################################################
S <- 100 ## spot
X <- 100 ## strike
tau <- 1 ## time-to-maturity
r <- 0.02 ## interest rate
q <- 0.08 ## dividend rate
v <- 0.2 ## volatility
## the closed-form solution
callBSM <- function(S,X,tau,r,q,v) {
d1 <- (log(S/X) + (r - q + v^2 / 2)*tau) / (v*sqrt(tau))
d2 <- d1 - v*sqrt(tau)
S * exp(-q * tau) * pnorm(d1) - X * exp(-r * tau) * pnorm(d2)
}
callBSM(S,X,tau,r,q,v)
## with the characteristic function
callCF(cf = cfBSM, S = S, X = X, tau = tau, r = r, q = q,
v = v^2, ## variance, not vol
implVol = TRUE)
###################################################
### code chunk number 128: NMOFman.Rnw:2601-2606
###################################################
cf <- c(5, 5, 5, 5, 5, 105) ## cashflows
times <- 1:6 ## times to payment
y <- 0.047 ## the "true" yield
b0 <- sum(cf/(1 + y)^times)
b0
###################################################
### code chunk number 129: NMOFman.Rnw:2615-2622
###################################################
vanillaBond <- function(cf, times, df, yields) {
if (missing(times))
times <- seq_len(length(cf))
if (missing(df))
df <- 1/(1+yields)^times
drop(cf %*% df)
}
###################################################
### code chunk number 130: NMOFman.Rnw:2627-2630
###################################################
cf <- c(rep(5, 9), 105)
vanillaBond(cf, yields = 0.05)
vanillaBond(cf, yields = 0.03)
###################################################
### code chunk number 131: NMOFman.Rnw:2639-2640
###################################################
2^(1:5)
###################################################
### code chunk number 132: NMOFman.Rnw:2651-2652
###################################################
vanillaBond(cf, 1:10, yield = NS(c(0.03,0,0,2), 1:10))
###################################################
### code chunk number 133: NMOFman.Rnw:2661-2667
###################################################
cf <- c(5, 5, 5, 5, 5, 105) ## cashflows
times <- 1:6 ## times to payment
y <- 0.047 ## the "true" yield
b0 <- sum(cf/(1 + y)^times)
cf <- c(-b0, cf); times <- c(0, times)
data.frame(times=times, cashflows=cf)
###################################################
### code chunk number 134: NMOFman.Rnw:2673-2709
###################################################
ytm <- function(cf, times, y0 = 0.05,
tol = 1e-05, h = 1e-05, maxit = 1000L) {
dr <- 1
for (i in seq_len(maxit)) {
y1 <- 1 + y0
g <- cf / y1 ^ times
g <- sum(g)
t1 <- times - 1
dg <- times * cf * 1/y1 ^ t1
dg <- sum(dg)
dr <- g/dg
y0 <- y0 + dr
if (abs(dr) < tol)
break
}
y0
}
ytm2 <- function(cf, times, y0 = 0.05,
tol = 1e-04, h = 1e-08, maxit = 1000L) {
dr <- 1
for (i in seq_len(maxit)) {
y1 <- 1 + y0
g <- sum(cf/y1^times)
y1 <- y1 + h
dg <- (sum(cf/y1^times) - g)/h
dr <- g/dg
y0 <- y0 - dr
if (abs(dr) < tol)
break
}
y0
}
system.time(for (i in 1:2000) ytm(cf, times, y0=0.06))
system.time(for (i in 1:2000) ytm2(cf, times, y0=0.06))
ytm(cf, times, y0=0.062, maxit = 5000)
ytm2(cf, times, y0=0.062, maxit = 5000)
###################################################
### code chunk number 135: NMOFman.Rnw:2720-2721
###################################################
(initial.value <- 5/b0)
###################################################
### code chunk number 136: NMOFman.Rnw:2724-2728
###################################################
ytm(cf, times, y0 = 0.7, maxit = 5000)
ytm(cf, times, y0 = initial.value)
ytm2(cf, times, y0 = 0.7, maxit = 5000)
ytm2(cf, times, y0 = initial.value)
###################################################
### code chunk number 137: NMOFman.Rnw:2774-2790
###################################################
mu <- 1
sigma <- 2
a <- -1
b <- 4
at <- (a - mu)/sigma
bt <- (b - mu)/sigma
## "throw away" strategy
x0 <- rnorm(10000L, mean = mu, sd = sigma)
x0 <- x0[x0>=a & x0<=b]
## inverse strategy
u <- runif(length(x0))
z <- qnorm(pnorm(at) + u*(pnorm(bt) - pnorm(at)))
x1 <- z * sigma + mu
###################################################
### code chunk number 138: NMOFman.Rnw:2795-2799
###################################################
par(mfrow = c(1, 2), mar = c(3, 3, 1, 1),
bty = "n", las = 1, ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
hist(x0, xlab = "")
###################################################
### code chunk number 139: NMOFman.Rnw:2801-2805
###################################################
par(mfrow = c(1, 2), mar = c(3, 3, 1, 1),
bty = "n", las = 1, ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
hist(x1, xlab = "")
###################################################
### code chunk number 140: NMOFman.Rnw:2811-2815
###################################################
x1 <- x1[1:750]
x2 <- rnorm(200)
x3 <- runif(500)
x4 <- rbinom(100, size = 50, prob = 0.4)
###################################################
### code chunk number 141: NMOFman.Rnw:2819-2821
###################################################
cormat <- array(0.5, dim = c(4, 4))
diag(cormat) <- 1
###################################################
### code chunk number 142: NMOFman.Rnw:2825-2828
###################################################
results <- resampleC(x1 = x1, x2 = x2, x3 = x3, x4 = x4,
size = 50, cormat = cormat)
cor(results, method = "spearman")
###################################################
### code chunk number 143: NMOFman.Rnw:2831-2845
###################################################
## this function is taken from ?pairs
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1L], y, col = grey(.5))
}
par(mar = c(3, 3, 1, 1),
bty = "n", las = 1, ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
pairs(results,
diag.panel = panel.hist,
gap = 0, pch = 19, cex = 0.5)
###################################################
### code chunk number 144: NMOFman.Rnw:2850-2860
###################################################
par(mfrow = c(2, 4), mar = c(3, 5, 1, 1),
bty = "n", las = 1, ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
hist(x1, xlab = "", ylab = "original")
hist(x2, xlab = "", ylab = "")
hist(x3, xlab = "", ylab = "")
hist(x4, xlab = "", ylab = "")
hist(results[ ,"x1"], xlab = "", ylab = "resampled")
hist(results[ ,"x2"], xlab = "", ylab = "")
hist(results[ ,"x3"], xlab = "", ylab = "")
hist(results[ ,"x4"], xlab = "", ylab = "")
###################################################
### code chunk number 145: NMOFman.Rnw:2943-2951
###################################################
na <- 500L ## number of assets
C <- array(0.6, dim = c(na,na)) ## correlation matrix
diag(C) <- 1
minVol <- 0.20; maxVol <- 0.40 ## covariance matrix
Vols <- (maxVol - minVol) * runif(na) + minVol
Sigma <- outer(Vols, Vols) * C
###################################################
### code chunk number 146: NMOFman.Rnw:2962-2968
###################################################
OF <- function(x, Data) {
sx <- sum(x)
w <- rep.int(1/sx, sx)
res <- crossprod(w, Data$Sigma[x, x])
tcrossprod(w, res)
}
###################################################
### code chunk number 147: NMOFman.Rnw:2971-2981
###################################################
neighbour <- function(xc, Data) {
xn <- xc
p <- sample.int(Data$na, Data$nn, replace = FALSE)
xn[p] <- !xn[p]
## reject infeasible solution
sumx <- sum(xn)
if ( (sumx > Data$Ksup) || (sumx < Data$Kinf) )
xc else xn
}
###################################################
### code chunk number 148: NMOFman.Rnw:2987-2992
###################################################
Data <- list(Sigma = Sigma, ## cov-matrix
Kinf = 30L, ## min cardinality
Ksup = 60L, ## max cardinality
na = na, ## number of assets
nn = 1L) ## how many assets to change per iteration
###################################################
### code chunk number 149: NMOFman.Rnw:2995-2999
###################################################
card0 <- sample(Data$Kinf:Data$Ksup, 1L, replace = FALSE)
assets <- sample.int(na, card0, replace = FALSE)
x0 <- logical(na)
x0[assets] <- TRUE
###################################################
### code chunk number 150: NMOFman.Rnw:3004-3012
###################################################
## Local Search
algo <- list(x0 = x0, neighbour = neighbour, nS = 5000L,
printDetail = FALSE, printBar = FALSE)
system.time(solLS <- LSopt(OF, algo = algo, Data = Data))
## Threshold Accepting
algo$nT <- 10L; algo$nS <- trunc(algo$nS/algo$nT); algo$q <- 0.2
system.time(solTA <- TAopt(OF, algo = algo, Data = Data))
###################################################
### code chunk number 151: NMOFman.Rnw:3096-3103
###################################################
OF2 <- function(x, Data) {
res <- colSums(Data$Sigma %*% x * x)
n <- colSums(x); res <- res / n^2
## penalise
p <- pmax(Data$Kinf - n, 0) + pmax(n - Data$Ksup, 0)
res + p
}
###################################################
### code chunk number 152: NMOFman.Rnw:3110-3113
###################################################
algo <- list(nB = na, nP = 100L, nG = 500L, prob = 0.002,
printBar = FALSE, loopOF = FALSE)
system.time(solGA <- GAopt(OF = OF2, algo = algo, Data = Data))
###################################################
### code chunk number 153: NMOFman.Rnw:3118-3122
###################################################
cat("Local Search ", format(sqrt(solLS$OFvalue), digits = 4), "\n",
"Threshold Accepting ", format(sqrt(solTA$OFvalue), digits = 4), "\n",
"Genetic Algorithm ", format(sqrt(solGA$OFvalue), digits = 4), "\n",
sep = "")
###################################################
### code chunk number 154: NMOFman.Rnw:3151-3159
###################################################
na <- 100L ## number of assets
ns <- 200L ## number of scenarios
vols <- runif(na, min = 0.2, max = 0.4) ## marginal vols
C <- matrix(0.6, na, na); diag(C) <- 1 ## correlation matrix
R <- rnorm(ns * na)/16 ## random returns
dim(R) <- c(ns, na)
R <- R %*% chol(C)
R <- R %*% diag(vols)
###################################################
### code chunk number 155: NMOFman.Rnw:3184-3194
###################################################
Data <- list(R = t(R), ## scenarios
theta = 0.005, ## return threshold
na = na, ## number of assets
ns = ns, ## number of scenarios
max = rep( 0.05, na), ## DE: vector of max. weight
min = rep(-0.05, na), ## DE: vector of min. weight
wsup = 0.05, ## TA: max weight
winf = -0.05, ## TA: min weight
eps = 0.5/100, ## TA: step size
w = 1) ## penalty weight
###################################################
### code chunk number 156: NMOFman.Rnw:3201-3204
###################################################
x0 <- Data$min + runif(Data$na)*(Data$max - Data$min)
x0[1:5]
sum(x0)
###################################################
### code chunk number 157: NMOFman.Rnw:3210-3214
###################################################
temp <- R %*% x0 ## compute portfolio returns
temp <- temp - Data$theta ## subtract return threshold
temp <- (temp[temp < 0])^2 ## select elements below threshold
sum(temp)/ns ## compute semivariance
###################################################
### code chunk number 158: NMOFman.Rnw:3220-3227
###################################################
OF <- function(x, Data) {
Rx <- crossprod(Data$R, x)
Rx <- Rx - Data$theta
Rx <- Rx - abs(Rx)
Rx <- Rx * Rx
colSums(Rx) /(4*Data$ns)
}
###################################################
### code chunk number 159: NMOFman.Rnw:3235-3237
###################################################
OF(x0, Data)
OF(cbind(x0, x0, x0), Data)
###################################################
### code chunk number 160: NMOFman.Rnw:3245-3258
###################################################
repair <- function(x, Data) {
myFun <- function(x)
x/sum(x)
if (is.null(dim(x)[2L]))
myFun(x) else apply(x, 2L, myFun)
}
repair2 <- function(x, Data) {
myFun <- function(x)
x + (1 - sum(x))/Data$na
if (is.null(dim(x)[2L]))
myFun(x) else apply(x, 2L, myFun)
}
###################################################
### code chunk number 161: NMOFman.Rnw:3264-3270
###################################################
sum(x0)
sum(repair(x0, Data))
sum(repair2(x0, Data))
colSums(repair( cbind(x0, x0, x0), Data))
colSums(repair2(cbind(x0, x0, x0), Data))
###################################################
### code chunk number 162: NMOFman.Rnw:3276-3278
###################################################
summary(repair (x0, Data)-x0)
summary(repair2(x0, Data)-x0)
###################################################
### code chunk number 163: NMOFman.Rnw:3283-3294
###################################################
penalty <- function(x, Data) {
up <- Data$max
lo <- Data$min
xadjU <- x - up
xadjU <- xadjU + abs(xadjU)
xadjL <- lo - x
xadjL <- xadjL + abs(xadjL)
if (is.null(dim(x)[2L]))
Data$w * (sum(xadjU) + sum(xadjL)) else
Data$w * (colSums(xadjU) + colSums(xadjL))
}
###################################################
### code chunk number 164: NMOFman.Rnw:3302-3308
###################################################
x0[1L] <- 0.30
penalty(x0, Data)
penalty(cbind(x0, x0, x0), Data)
x0[1L] <- 0
penalty(x0, Data)
penalty(cbind(x0, x0, x0), Data)
###################################################
### code chunk number 165: NMOFman.Rnw:3313-3326
###################################################
algo <- list(nP = 100, ## population size
nG = 1000, ## number of generations
F = 0.25, ## step size
CR = 0.9,
min = Data$min,
max = Data$max,
repair = repair,
pen = penalty,
printBar = FALSE,
printDetail = TRUE,
loopOF = TRUE, ## do not vectorise
loopPen = TRUE, ## do not vectorise
loopRepair = TRUE) ## do not vectorise
###################################################
### code chunk number 166: NMOFman.Rnw:3332-3339
###################################################
system.time(sol <- DEopt(OF = OF,algo = algo,Data = Data))
16 * 100 * sqrt(sol$OFvalue) ## solution quality
## check constraints
all(all.equal(sum(sol$xbest), 1), ## budget constraint
sol$xbest <= Data$max, ## holding size constraints
sol$xbest >= Data$min)
###################################################
### code chunk number 167: NMOFman.Rnw:3349-3359
###################################################
## looping over the population
algo$loopOF <- TRUE; algo$loopPen <- TRUE; algo$loopRepair <- TRUE
t1 <- system.time(sol <- DEopt(OF = OF,algo = algo, Data = Data))
## evaluating the population in one step
algo$loopOF <- FALSE; algo$loopPen <- FALSE; algo$loopRepair <- FALSE
t2 <- system.time(sol <- DEopt(OF = OF,algo = algo, Data = Data))
## speedup
t1[[3L]]/t2[[3L]]
###################################################
### code chunk number 168: NMOFman.Rnw:3370-3383
###################################################
algo$printDetail <- FALSE
restartsDE <- restartOpt(fun = DEopt, ## what function
n = 20L, ## how many restarts
OF = OF,
algo = algo,
Data = Data,
method = "snow", ## using package snow
cl = 2) ## 2 cores
## extract best solution
OFvaluesDE <- sapply(restartsDE, `[[`, "OFvalue")
OFvaluesDE <- 16 * 100 * sqrt(OFvaluesDE)
weightsDE <- sapply(restartsDE, `[[`, "xbest")
###################################################
### code chunk number 169: NMOFman.Rnw:3388-3393
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0),
ps = 8, tck = 0.001)
plot(sort(OFvaluesDE), (seq_len(length(OFvaluesDE))) / length(OFvaluesDE),
type = "S", ylim = c(0, 1), xlab = "", ylab = "")
mtext("OF value", side = 1, line = 2)
###################################################
### code chunk number 170: NMOFman.Rnw:3397-3403
###################################################
par(bty = "n", las = 1, mar = c(3, 4, 0, 0),
ps = 8, tck = 0.001)
boxplot(t(weightsDE),
outline = FALSE, boxwex = 0.4, ylim = c(-0.06,0.06))
mtext("assets", side = 1, line = 2)
mtext("weights", side = 2, line = 1.3, las = 1, padj = -4)
###################################################
### code chunk number 171: NMOFman.Rnw:3416-3431
###################################################
algo$printDetail <- FALSE; algo$nP <- 200L; restarts <- 20L
nGs <- c(500L, 1500L, 3000L)
lstOFvaluesDE <- list()
for (i in 1:3) {
algo$nG <- nGs[i]
restartsDE <- restartOpt(fun = DEopt,
n = restarts,
OF = OF, algo = algo, Data = Data,
method = "snow", cl = 2)
## extract best solution
OFvaluesDE <- sapply(restartsDE, `[[`, "OFvalue")
OFvaluesDE <- 16 * 100 * sqrt(OFvaluesDE)
lstOFvaluesDE[[i]] <- OFvaluesDE
}
res <- simplify2array(lstOFvaluesDE)
###################################################
### code chunk number 172: NMOFman.Rnw:3436-3450
###################################################
algo$repair <- repair2
lstOFvaluesDE <- list()
for (i in 1:3) {
algo$nG <- nGs[i]
restartsDE <- restartOpt(fun = DEopt,
n = restarts,
OF = OF, algo = algo, Data = Data,
method = "snow", cl = 2)
## extract best solution
OFvaluesDE <- sapply(restartsDE, `[[`, "OFvalue")
OFvaluesDE <- 16 * 100 * sqrt(OFvaluesDE)
lstOFvaluesDE[[i]] <- OFvaluesDE
}
res2 <- simplify2array(lstOFvaluesDE)
###################################################
### code chunk number 173: NMOFman.Rnw:3455-3465
###################################################
allres <- as.vector(rbind(res,res2))
xlims <- pretty(allres); xlims <- c(min(xlims), max(xlims))
par(bty = "n", las = 1, mar = c(3, 4, 0, 0),
ps = 8, tck = 0.001)
plot(ecdf(res[ ,3L]), xlim = xlims, cex = 0.4,
main = "", ylab = "", xlab = "")
for (i in 1:2)
lines(ecdf(res[ ,i]), cex = 0.4)
for (i in 1:3)
lines(ecdf(res2[ ,i]), col = "blue", cex = 0.4)
###################################################
### code chunk number 174: NMOFman.Rnw:3475-3482
###################################################
weightsDE <- sapply(restartsDE, `[[`, "xbest")
par(bty = "n", las = 1, mar = c(3, 4, 0, 0),
ps = 8, tck = 0.001)
boxplot(t(weightsDE),
outline = FALSE, boxwex = 0.4, ylim = c(-0.06, 0.06))
mtext("assets", side = 1, line = 2)
mtext("weights", side = 2, line = 1.3, las = 1, padj = -4)
###################################################
### code chunk number 175: NMOFman.Rnw:3496-3513
###################################################
algo <- list(nP = 100L, ## population size
nG = 1000L, ## number of generations
c1 = 0.5, ## weight for individually best solution
c2 = 1.5, ## weight for overall best solution
min = Data$min,
max = Data$max,
repair = repair, pen = penalty,
iner = 0.7, initV = 1, maxV = 0.2,
printBar = FALSE, printDetail = TRUE)
system.time(sol <- PSopt(OF = OF,algo = algo,Data = Data))
16 * 100 * sqrt(sol$OFvalue) ## solution quality
## check constraints
all(all.equal(sum(sol$xbest),1), ## budget constraint
sol$xbest <= Data$max,
sol$xbest >= Data$min)
###################################################
### code chunk number 176: NMOFman.Rnw:3521-3528
###################################################
changeV <- function(x, Data) {
myFun <- function(x) x - (sum(x))/Data$na
if (is.null(dim(x)[2L]))
myFun(x) else apply(x, 2L, myFun)
}
sum(changeV(x0, Data))
colSums(changeV(cbind(x0, x0, x0), Data))
###################################################
### code chunk number 177: NMOFman.Rnw:3533-3537
###################################################
initP <- Data$min + diag(Data$max - Data$min) %*%
array(runif(length(Data$min) * algo$nP),
dim = c(length(Data$min), algo$nP))
colSums(initP <- repair(initP,Data))[1:10] ## check
###################################################
### code chunk number 178: NMOFman.Rnw:3542-3548
###################################################
algo$changeV <- changeV ## function to adjust velocity
algo$initP <- initP ## initial population
algo$repair <- NULL ## not needed anymore
system.time(sol <- PSopt(OF = OF,algo = algo, Data = Data))
16 * 100 * sqrt(sol$OFvalue) ## solution quality
###################################################
### code chunk number 179: NMOFman.Rnw:3552-3555
###################################################
all(all.equal(sum(sol$xbest), 1), ## budget constraint
sol$xbest <= Data$max,
sol$xbest >= Data$min)
###################################################
### code chunk number 180: NMOFman.Rnw:3558-3561
###################################################
algo$loopOF <- FALSE; algo$loopPen <- FALSE
algo$loopRepair <- FALSE; algo$loopChangeV <- FALSE
system.time(sol <- PSopt(OF = OF, algo = algo, Data = Data))
###################################################
### code chunk number 181: NMOFman.Rnw:3566-3582
###################################################
algo$printDetail <- FALSE
restartsPS <- restartOpt(fun = PSopt,
n = 20L,
OF = OF,
algo = algo, Data = Data,
method = "snow", cl = 2)
## extract best solution
OFvaluesPS <- sapply(restartsPS, `[[`, "OFvalue")
OFvaluesPS <- 16 * 100 * sqrt(OFvaluesPS)
par(bty = "n", las = 1,mar = c(3,4,0,0),
ps = 8, tck = 0.001)
plot(sort(OFvaluesPS),
(seq_len(length(OFvaluesPS))) / length(OFvaluesPS),
type = "S", ylim = c(0, 1), xlab = "", ylab = "")
mtext("OF value", side = 1, line = 2)
###################################################
### code chunk number 182: NMOFman.Rnw:3593-3615
###################################################
Data$R <- R ## not transposed any more
neighbourU <- function(sol, Data){
resample <- function(x, ...)
x[sample.int(length(x), ...)]
wn <- sol$w
toSell <- wn > Data$winf
toBuy <- wn < Data$wsup
i <- resample(which(toSell), size = 1L)
j <- resample(which(toBuy), size = 1L)
eps <- runif(1) * Data$eps
eps <- min(wn[i] - Data$winf, Data$wsup - wn[j], eps)
wn[i] <- wn[i] - eps
wn[j] <- wn[j] + eps
Rw <- sol$Rw + Data$R[,c(i,j)] %*% c(-eps,eps)
list(w = wn, Rw = Rw)
}
OF <- function(x, Data) {
Rw <- x$Rw - Data$theta
Rw <- Rw - abs(Rw)
sum(Rw*Rw) / (4*Data$ns)
}
###################################################
### code chunk number 183: NMOFman.Rnw:3619-3631
###################################################
w0 <- runif(Data$na); w0 <- w0/sum(w0)
x0 <- list(w = w0, Rw = R %*% w0)
algo <- list(x0 = x0,
neighbour = neighbourU,
nS = 2000L,
nT = 10L,
nD = 5000L,
q = 0.20,
printBar = FALSE,
printDetail = FALSE)
system.time(sol2 <- TAopt(OF,algo,Data))
16 * 100 * sqrt(sol2$OFvalue)
###################################################
### code chunk number 184: NMOFman.Rnw:3635-3657
###################################################
restartsTA <- restartOpt(fun = TAopt,
n = 20L,
OF = OF,
algo = algo,
Data = Data,
method = "snow",
cl = 2)
OFvaluesTA <- sapply(restartsTA, `[[`, "OFvalue") ## extract best solution
OFvaluesTA <- 16 * 100 * sqrt(OFvaluesTA)
weightsTA <- sapply(restartsTA, `[[`, "xbest")
par(bty = "n", las = 1,mar = c(3,4,0,0), ps = 8,
tck = 0.001, mgp = c(3, 0.5, 0))
## blue: DE solution with nP = 200 and nG = 2000
xlims <- pretty(c(res2[,3], OFvaluesTA))
plot(ecdf(res2[,3]), col = "blue", cex = 0.4,
main = "", ylab = "", xlab = "",
xlim = c(min(xlims), max(xlims)) )
## black: TA
lines(ecdf(OFvaluesTA), cex = 0.4)
###################################################
### code chunk number 185: NMOFman.Rnw:3680-3686
###################################################
na <- 50
no <- 5000
D <- array(rnorm(na*no)*0.01, dim = c(no,na))
w <- runif(na)
w <- w/sum(w)
R <- D %*% w
###################################################
### code chunk number 186: NMOFman.Rnw:3691-3692
###################################################
require("compiler")
###################################################
### code chunk number 187: NMOFman.Rnw:3699-3705
###################################################
var1 <- function(R) {
n <- NROW(R)
m <- sum(R)/n
crossprod(R)/(n-1) - m^2
}
var(R) - var1(R)
###################################################
### code chunk number 188: NMOFman.Rnw:3711-3724
###################################################
var1 <- function(R) {
n <- NROW(R)
m <- sum(R)/n
crossprod(R)/(n - 1) - m^2
}
var2 <- cmpfun(var1)
runs <- 10000
system.time(for (i in seq_len(runs))
ignore <- var(R))
system.time(for (i in seq_len(runs))
ignore <- var1(R))
system.time(for (i in seq_len(runs))
ignore <- var2(R))
###################################################
### code chunk number 189: NMOFman.Rnw:3731-3758
###################################################
pm0 <- function(x, xp = 2, threshold = 0, lower = TRUE) {
n <- NROW(x)
x <- x - threshold
if (lower)
x <- x[x < 0] else x <- x[x > 0]
sum(x^xp)/n
}
pm1 <- function(x, xp = 2, threshold = 0, lower = TRUE, keep.sign = FALSE) {
x <- x - threshold
if (lower)
x <- x - abs(x)
else
x <- x + abs(x)
sx <- sign(x)
x <- abs(x)
if (xp == 1L)
sum(x)/2/length(x)
else if (xp == 2L)
sum(x*x)/4/length(x)
else if (xp == 3L)
sum(x*x*x)/8/length(x)
else if (xp == 4L)
sum(x*x*x*x)/16/length(x)
else
sum(x^xp)/2^xp/length(x)
}
pm2 <- cmpfun(pm1)
###################################################
### code chunk number 190: NMOFman.Rnw:3764-3785
###################################################
pm0(R)
pm1(R)
pm2(R)
runs <- 1000
system.time(for (i in seq_len(runs))
ignore <- pm0(R))
system.time(for (i in seq_len(runs))
ignore <- pm1(R))
system.time(for (i in seq_len(runs))
ignore <- pm2(R))
pm0(R,2.2)
pm1(R,2.2)
pm2(R,2.2)
runs <- 1000
system.time(for (i in seq_len(runs))
ignore <- pm0(R, 2.5))
system.time(for (i in seq_len(runs))
ignore <- pm1(R, 2.5))
system.time(for (i in seq_len(runs))
ignore <- pm2(R, 2.5))
###################################################
### code chunk number 191: NMOFman.Rnw:3822-3829
###################################################
size <- 20L
x <- logical(size)
x[runif(size) > 0.5] <- TRUE
## store information
Data <- list()
Data$size <- size
###################################################
### code chunk number 192: NMOFman.Rnw:3833-3853
###################################################
compareLogicals <- function(x,y, sep = "", ## true = "1", false = "0",
mark = "^", below = TRUE) {
mark.line <- ifelse(x == y, " ", mark)
if (!below)
cat(mark.line, "\n", sep = sep)
cat(as.integer(x), "\n",
as.integer(y), "\n", sep = sep)
if (below)
cat(mark.line, "\n", sep = sep)
sxy <- sum(x != y)
if (!sxy)
cat("The vectors do not differ.\n", sep = "")
else if (sxy == 1L)
cat("The vectors differ in 1 place.\n", sep = "")
else
cat("The vectors differ in ", sum(x != y), " place(s).\n", sep = "")
invisible(x != y)
}
###################################################
### code chunk number 193: NMOFman.Rnw:3860-3861
###################################################
compareLogicals(x, x) ## there should be no difference
###################################################
### code chunk number 194: NMOFman.Rnw:3865-3867
###################################################
z <- x; z[2L] <- !z[2L]
compareLogicals(x, z)
###################################################
### code chunk number 195: NMOFman.Rnw:3876-3883
###################################################
Data$n <- 5L ## how many elements to change
neighbour <- function(x, Data) {
ii <- sample.int(Data$size, Data$n)
x[ii] <- !x[ii]
x
}
compareLogicals(x, neighbour(x, Data))
###################################################
### code chunk number 196: NMOFman.Rnw:3893-3904
###################################################
neighbour <- function(x, Data) {
## required: x must have at least one TRUE and one FALSE
Ts <- which(x)
Fs <- which(!x)
lenTs <- length(Ts)
O <- sample.int(lenTs, 1L)
I <- sample.int(Data$size - lenTs, 1L)
x[c(Fs[I], Ts[O])] <- c(TRUE, FALSE)
x
}
compareLogicals(x, neighbour(x, Data))
###################################################
### code chunk number 197: NMOFman.Rnw:3911-3929
###################################################
size <- 5L
x0 <- runif(size)
xTRUE <- runif(size)
Data <- list(xTRUE = xTRUE,
step = 0.02)
OF <- function(x, Data)
max(abs(x - Data$xTRUE))
neighbour <- function(x, Data)
x + runif(length(Data$xTRUE))*Data$step - Data$step/2
algo <- list(q = 0.05, nS = 1000L, nT = 10L,
neighbour = neighbour, x0 = x0,
printBar = FALSE,
printDetail = FALSE,
storeSolutions = TRUE,
storeF = TRUE)
res <- TAopt(OF, algo = algo, Data = Data)
res$OFvalue < 0.005
###################################################
### code chunk number 198: NMOFman.Rnw:3950-4221
###################################################
## N1: This neighbour enforces a budget constraint, a minimum
## holding size (which need not be zero) and a maximum holding size
Data <- list(wmin = 0, ## the minimal weight
wmax = 0.22, ## the maximal weight
eps = 0.2/100, ## the step size
## resample = function(x, ...)
## x[sample.int(length(x), ...)],
na = dim(fundData)[2L],
R = fundData)
cat("The portfolio will consist of at least ",
ceiling(1/Data$wmax), " assets.\n", sep = "")
neighbour1 <- function(w, Data){
toSell <- which(w > Data$wmin)
toBuy <- which(w < Data$wmax)
i <- toSell[sample.int(length(toSell), size = 1L)]
j <- toBuy[ sample.int(length(toBuy), size = 1L)]
eps <- runif(1) * Data$eps
eps <- min(w[i] - Data$wmin, Data$wmax - w[j], eps)
w[i] <- w[i] - eps
w[j] <- w[j] + eps
w
}
neighbour1U <- function(x, Data){
wn <- x$w
toSell <- which(wn > Data$wmin)
toBuy <- which(wn < Data$wmax)
i <- toSell[sample.int(length(toSell), size = 1L)]
j <- toBuy[ sample.int(length(toBuy), size = 1L)]
eps <- runif(1) * Data$eps
eps <- min(wn[i] - Data$wmin, Data$wmax - wn[j], eps)
wn[i] <- wn[i] - eps
wn[j] <- wn[j] + eps
Rw <- x$Rw + Data$R[ ,c(i,j)] %*% c(-eps,eps)
list(w = wn, Rw = Rw)
}
## create a random solution
makex <- function(Data) {
resample <- function(x, ...)
x[sample.int(length(x), ...)]
w0 <- numeric(Data$na)
nAssets <- resample(ceiling(1/Data$wmax):Data$na, 1L)
w0[sample(seq_len(Data$na), nAssets)] <- runif(nAssets)
w0/sum(w0)
}
isOK <- function(w, Data) {
tooBig <- any(w > Data$wmax)
tooSmall <- any(w < Data$wmin)
sumToOne <- abs(sum(w)-1) < 1e-12
if (!tooBig && !tooSmall && sumToOne)
TRUE
else
FALSE
}
## TEST 1
w0 <- makex(Data)
x0 <- list(w = w0, Rw = fundData %*% w0)
isOK(w0, Data)
isOK(x0$w, Data)
set.seed(545)
w0 <- makex(Data)
nTests <- 1e3
for (i in seq(nTests)) {
w1 <- neighbour1(w0, Data)
if (isOK(w1, Data))
w0 <- w1
else
stop("error")
}
set.seed(545)
w0 <- makex(Data)
x0 <- list(w = w0, Rw = fundData %*% w0)
nTests <- 1e3
for (i in seq(nTests)) {
x1 <- neighbour1U(x0, Data)
if (isOK(x1$w, Data))
x0 <- x1
else
stop("error")
}
all.equal(fundData %*% w1, x1$Rw)
## TEST 2: reach a target solution
makeOF <- function(wt)
function(w0, Data)
sum(abs(wt - w0))
wt <- makex(Data)
OF <- makeOF(wt)
w0 <- makex(Data)
OF(w0, Data)
TAsettings <- list(neighbour = neighbour1,
x0 = w0, nS = 5000, q = 0.1,
printBar = FALSE)
res <- TAopt(OF, algo = TAsettings, Data)
round(head(sort(abs(res$xbest-wt), decreasing = TRUE),5),6)
## N2: This long-only neighbour enforces a budget constraint, a
## minimum holding size (which need not be zero), and a maximum holding
## size and a maximum cardinality.
Data <- list(wmax = 0.3, ## the maximal weight
Kmax = 10, ## max cardinality
eps = 1/100, ## the step size
## resample = function(x, ...)
## x[sample.int(length(x), ...)],
na = dim(fundData)[2L],
R = fundData)
cat("The portfolio will consist of at least ",
ceiling(1/Data$wmax), " assets.\n", sep = "")
neighbour2 <- function(w, Data){
tol <- 1e-12
J <- sum(w > tol)
if (J == Data$Kmax)
toBuy <- which(w > tol & w < Data$wmax)
else
toBuy <- which(w < Data$wmax)
toSell <- which(w > tol)
i <- toSell[sample.int(length(toSell), size = 1L)]
j <- toBuy[ sample.int(length(toBuy), size = 1L)]
eps <- runif(1) * Data$eps
eps <- min(w[i], Data$wmax - w[j], eps)
w[i] <- w[i] - eps
w[j] <- w[j] + eps
w
}
neighbour2U <- function(x, Data){
tol <- 1e-12
w <- x$w
J <- sum(w > tol)
if (J == Data$Kmax)
toBuy <- which(w > tol & w < Data$wmax)
else
toBuy <- which(w < Data$wmax)
toSell <- which(w > tol)
i <- toSell[sample.int(length(toSell), size = 1L)]
j <- toBuy[ sample.int(length(toBuy), size = 1L)]
eps <- runif(1) * Data$eps
eps <- min(w[i], Data$wmax - w[j], eps)
w[i] <- w[i] - eps
w[j] <- w[j] + eps
Rw <- x$Rw + Data$R[ ,c(i,j)] %*% c(-eps, eps)
list(w = w, Rw = Rw)
}
makex <- function(Data) {
w0 <- numeric(Data$na)
nAssets <- sample(ceiling(1/Data$wmax):Data$Kmax, 1L)
w0[sample(seq_len(Data$na), nAssets)] <- runif(nAssets)
w0/sum(w0)
}
isOK <- function(w, Data) {
tooBig <- any(w > Data$wmax)
tooMany <- sum(w > 1e-12) > Data$Kmax
sumToOne <- abs(sum(w)-1) < 1e-12
if (!tooBig && !tooMany && sumToOne)
TRUE
else
FALSE
}
## TEST 1
w0 <- makex(Data)
x0 <- list(w = w0, Rw = fundData %*% w0)
isOK(w0, Data)
isOK(x0$w, Data)
set.seed(545)
w0 <- makex(Data)
nTests <- 1e3
for (i in seq(nTests)) {
w1 <- neighbour2(w0, Data)
if (isOK(w1, Data))
w0 <- w1
else
stop("error")
}
set.seed(545)
w0 <- makex(Data)
x0 <- list(w = w0, Rw = fundData %*% w0)
nTests <- 1e3
for (i in seq(nTests)) {
x1 <- neighbour2U(x0, Data)
if (isOK(x1$w, Data))
x0 <- x1
else
stop("error")
}
all.equal(fundData %*% w1, x1$Rw)
## TEST 2: reach a target solution
makeOF <- function(wt)
function(w0, Data)
sum(abs(wt - w0))
wt <- makex(Data)
OF <- makeOF(wt)
w0 <- makex(Data)
OF(w0, Data)
OF(wt, Data)
TAsettings <- list(neighbour = neighbour2,
x0 = w0, nS = 5000, q = 0.1,
printBar = FALSE)
res <- TAopt(OF, algo = TAsettings, Data)
isOK(res$xbest, Data)
df <- data.frame(target=wt, w0 = w0, wTAopt = res$xbest)
tmpfun <- function(x)
!all(abs(x) < 1e-14)
df[apply(df,1,tmpfun),]
apply(df, 2, sum)
wt <- numeric(200)
wt[1:4] <- c(0.3,0.3,0.3,0.1)
OF <- makeOF(wt)
TAsettings <- list(neighbour = neighbour2,
x0 = w0, nS = 5000, q = 0.1,
printBar = FALSE)
res <- TAopt(OF, algo = TAsettings, Data)
isOK(res$xbest, Data)
df <- data.frame(target=wt, w0 = w0, wTAopt = res$xbest)
tmpfun <- function(x)
!all(abs(x) < 1e-14)
df[apply(df,1,tmpfun),]
apply(df, 2, sum)
w0 <- makex(Data)
x0 <- list(w = w0, Rw = fundData %*% w0)
## the N is slower
system.time(for (i in 1:10000) neighbour2(w0, Data))
system.time(for (i in 1:10000) neighbour2U(x0, Data))
TAsettings2 <- list(neighbour = neighbour2,
x0 = w0, nS = 500, q = 0.1,
printBar = FALSE, printDetail = FALSE)
TAsettings2U <- list(neighbour = neighbour2U,
x0 = x0, nS = 500, q = 0.1,
printBar = FALSE, printDetail = FALSE)
ofun <- function(w, Data) {
Rw <- Data$R %*% w
crossprod(Rw)
}
ofunU <- function(sol, Data)
crossprod(sol$Rw)
ign <- TAopt(ofun, TAsettings2, Data)
ign <- TAopt(ofunU, TAsettings2U, Data)
##benchmark(ign <- TAopt(ofun, TAsettings2, Data),
## ign <- TAopt(ofunU, TAsettings2U, Data),
## replications = 1, order = "relative")
###################################################
### code chunk number 199: NMOFman.Rnw:4256-4258
###################################################
maxabs <- function(x, lim)
max(sum(abs(x)) - lim, 0)
###################################################
### code chunk number 200: NMOFman.Rnw:4263-4346
###################################################
require("compiler")
wmin <- 0.01
wmax <- 0.10
w <- numeric(50)
w[1:10] <- 0.1
step <- 0.01
N <- function(w) {
wo <- w
## initial sell
sel <- which(w >= wmin)
i <- sel[sample.int(length(sel), size = 1)]
if (w[i] == wmin) {
eps <- wmin
w[i] <- 0
} else {
eps <- min(runif(1)*step , w[i] - wmin)
w[i] <- w[i] - eps
}
cash <- eps
iter <- 0
while (abs(cash) > 1e-14) {
iter <- iter + 1
if (iter > 10) {
return(wo)
}
##message("cash ", cash)
if (cash > 0) { ## buy something
sel <- which(w < wmax)
i <- sel[sample.int(length(sel), size = 1)]
if (w[i] == 0) {
w[i] <- eps <- wmin
} else {
eps <- min(runif(1)*step , wmax - w[i], cash)
w[i] <- w[i] + eps
}
cash <- cash - eps
} else { ## sell something
sel <- which(w >= wmin)
i <- sel[sample.int(length(sel), size = 1)]
if (w[i] == wmin) {
eps <- wmin
w[i] <- 0
} else {
eps <- min(runif(1)*step , w[i] - wmin)
w[i] <- w[i] - eps
}
cash <- cash + eps
}
}
##message(iter)
w
}
##N <- cmpfun(N)
##w
system.time(for (i in 1:10000) w <- N(w))
## goal <- numeric(50)
## goal[41:50] <- 0.1
## OF <- function(x) {
## tmp <- x - goal
## sum(tmp * tmp)
## }
## ans <- LSopt(OF, algo = list(nS = 1000000, neighbour = N, x0 = w))
## sum(ans$xbest > 0)
## ans$xbest
## ans <- numeric(10000)
## for (i in seq_along(ans))
## ans[i] <- ceiling(runif(1)*10)
## system.time(for (i in 1:10000) ignore <- ceiling(runif(5)*99))
## system.time(for (i in 1:10000) ignore <- sample.int(99, 5))
###################################################
### code chunk number 201: NMOFman.Rnw:4361-4372
###################################################
makeCashFlows <- function(coupon, T) {
t1 <- T - floor(T) ## time to first coupon
tm <- seq(ifelse(t1 > 1e-5, t1, 1), ## 1e-5 is less than a calendar day
T,
by = 1)
cf <- rep.int(coupon, length(tm))
cf[length(cf)] <- cf[length(cf)] + 100
list(cf = cf, tm = tm)
}
makeCashFlows(3, 10.2)
###################################################
### code chunk number 202: NMOFman.Rnw:4376-4396
###################################################
cf1 <- c(rep(5.75, 8), 105.75); tm1 <- 0:8 + 0.5
cf2 <- c(rep(4.25, 17), 104.25); tm2 <- 1:18
cf3 <- c(3.5, 103.5); tm3 <- 0:1 + 0.5
cf4 <- c(rep(3.00, 15), 103.00); tm4 <- 1:16
cf5 <- c(rep(3.25, 11), 103.25); tm5 <- 0:11 + 0.5
cf6 <- c(rep(5.75, 17), 105.75); tm6 <- 0:17 + 0.5
cf7 <- c(rep(3.50, 14), 103.50); tm7 <- 1:15
cf8 <- c(rep(5.00, 8), 105.00); tm8 <- 0:8 + 0.5
cf9 <- 105; tm9 <- 1
cf10 <- c(rep(3.00, 12), 103.00); tm10 <- 0:12 + 0.5
cf11 <- c(rep(2.50, 7), 102.50); tm11 <- 1:8
cf12 <- c(rep(4.00, 10), 104.00); tm12 <- 1:11
cf13 <- c(rep(3.75, 18), 103.75); tm13 <- 0:18 + 0.5
cf14 <- c(rep(4.00, 17), 104.00); tm14 <- 1:18
cf15 <- c(rep(2.25, 8), 102.25); tm15 <- 0:8 + 0.5
cf16 <- c(rep(4.00, 6), 104.00); tm16 <- 1:7
cf17 <- c(rep(2.25, 12), 102.25); tm17 <- 1:13
cf18 <- c(rep(4.50, 19), 104.50); tm18 <- 0:19 + 0.5
cf19 <- c(rep(2.25, 7), 102.25); tm19 <- 1:8
cf20 <- c(rep(3.00, 14), 103.00); tm20 <- 1:15
###################################################
### code chunk number 203: NMOFman.Rnw:4402-4417
###################################################
cfList <- list( cf1, cf2, cf3, cf4, cf5, cf6, cf7, cf8, cf9,cf10,
cf11,cf12,cf13,cf14,cf15,cf16,cf17,cf18,cf19,cf20)
tmList <- list( tm1, tm2, tm3, tm4, tm5, tm6, tm7, tm8, tm9,tm10,
tm11,tm12,tm13,tm14,tm15,tm16,tm17,tm18,tm19,tm20)
tm <- unlist(tmList, use.names = FALSE)
tm <- sort(unique(tm))
nR <- length(tm)
nC <- length(cfList)
cfMatrix <- array(0, dim = c(nR, nC))
for(j in seq(nC))
cfMatrix[tm %in% tmList[[j]], j] <- cfList[[j]]
rownames(cfMatrix) <- tm
cfMatrix[1:10, 1:10]
###################################################
### code chunk number 204: NMOFman.Rnw:4424-4428
###################################################
betaTRUE <- c(5,-2,1,10,1,3)
yM <- NSS(betaTRUE,tm)
diFa <- 1 / ( (1 + yM/100)^tm )
bM <- diFa %*% cfMatrix
###################################################
### code chunk number 205: NMOFman.Rnw:4432-4436
###################################################
Data <- list(bM = bM, tm = tm, cfMatrix = cfMatrix, model = NSS,
ww = 1,
min = c( 0,-15,-30,-30,0 ,2.5),
max = c(15, 30, 30, 30,2.5,5 ))
###################################################
### code chunk number 206: NMOFman.Rnw:4446-4456
###################################################
OF2 <- function(param, Data) {
tm <- Data$tm
bM <- Data$bM
cfMatrix <- Data$cfMatrix
diFa <- 1 / ((1 + Data$model(param, tm)/100)^tm)
b <- diFa %*% cfMatrix
aux <- b - bM; aux <- max(abs(aux))
if (is.na(aux)) aux <- 1e10
aux
}
###################################################
### code chunk number 207: NMOFman.Rnw:4460-4475
###################################################
penalty <- function(mP, Data) {
minV <- Data$min
maxV <- Data$max
ww <- Data$ww
## if larger than maxV, element in A is positiv
A <- mP - as.vector(maxV)
A <- A + abs(A)
## if smaller than minV, element in B is positiv
B <- as.vector(minV) - mP
B <- B + abs(B)
## beta 1 + beta2 > 0
C <- ww*((mP[1L, ] + mP[2L, ]) - abs(mP[1L, ] + mP[2L, ]))
A <- ww * colSums(A + B) - C
A
}
###################################################
### code chunk number 208: NMOFman.Rnw:4480-4496
###################################################
algo <- list(nP = 200L,
nG = 1000L,
F = 0.50,
CR = 0.99,
min = c( 0,-15,-30,-30,0 ,2.5),
max = c(15, 30, 30, 30,2.5,5 ),
pen = penalty,
repair = NULL,
loopOF = TRUE,
loopPen = FALSE,
loopRepair = FALSE,
printBar = FALSE,
printDetail = FALSE,
storeF = FALSE)
sol <- DEopt(OF = OF2, algo = algo, Data = Data)
###################################################
### code chunk number 209: NMOFman.Rnw:4503-4505
###################################################
max( abs(Data$model(sol$xbest, tm) - Data$model(betaTRUE, tm)))
sol$OFvalue
###################################################
### code chunk number 210: NMOFman.Rnw:4510-4527
###################################################
s0 <- algo$min + (algo$max - algo$min) * runif(length(algo$min))
system.time(sol2 <- nlminb(s0,OF2,Data = Data,
lower = Data$min,
upper = Data$max,
control = list(eval.max = 50000,
iter.max = 50000)))
max(abs(Data$model(sol2$par,tm) - Data$model(betaTRUE,tm)))
sol2$objective
par(ps = 8, bty = "n", las = 1, tck = 0.01,
mgp = c(3, 0.5, 0), mar = c(4, 4, 1, 1))
plot(tm, yM, xlab = "maturities in years", ylab = "yields in %")
lines(tm,Data$model(sol$xbest,tm), col = "blue")
lines(tm,Data$model(sol2$par,tm), col = "darkgreen", lty = 2)
legend(x = "bottom", legend = c("true yields", "DE", "nlminb"),
col = c("black", "blue", "darkgreen"),
pch = c(1, NA, NA), lty = c(0, 1, 2))
###################################################
### code chunk number 211: NMOFman.Rnw:4532-4535
###################################################
diFa <- 1 / ((1 + NSS(sol$xbest,tm)/100)^tm)
b <- diFa %*% cfMatrix
b - bM
###################################################
### code chunk number 212: NMOFman.Rnw:4540-4544
###################################################
par(ps = 8, bty = "n", las = 1, tck = 0.01,
mgp = c(3, 0.5, 0), mar = c(4, 4, 1, 1))
plot(tm, NSS(sol$xbest,tm) - NSS(betaTRUE,tm),
xlab = "maturities in years", ylab = "yield error in %")
###################################################
### code chunk number 213: NMOFman.Rnw:4551-4555
###################################################
par(ps = 8, bty = "n", las = 1, tck = 0.01,
mgp = c(3, 0.5, 0), mar = c(4, 4, 1, 1))
plot(as.numeric(unlist(lapply(tmList, max))), as.vector(b - bM),
xlab = "maturities in years", ylab = "price error in %")
###################################################
### code chunk number 214: NMOFman.Rnw:4563-4567
###################################################
beta <- c(5,-2,1,10,1,3)
yM <- NSS(beta,tm)
diFa <- 1 / ( (1 + yM/100)^tm )
b <- diFa %*% cfMatrix
###################################################
### code chunk number 215: NMOFman.Rnw:4572-4578
###################################################
B <- cbind(c(5,-2,1,10,1,3), c(4,-2,1,10,1,3))
Y <- array(0, dim = c(length(tm), ncol(B)))
for (i in 1:ncol(Y))
Y[ ,i] <- NSS(B[ ,i], tm)
D <- 1/((1+Y/100)^tm)
t(cfMatrix) %*% D - as.vector(b)
###################################################
### code chunk number 216: NMOFman.Rnw:4589-4608
###################################################
fy <- function(ytm, cf, tm)
sum( cf / ( (1 + ytm)^tm ) )
compYield <- function(cf, tm, guess = NULL) {
logik <- cf != 0
cf <- cf[logik]
tm <- tm[logik]
if (is.null(guess)) {ytm <- 0.05} else {ytm <- guess}
h <- 1e-8; dF <- 1; ci <- 0
while (abs(dF) > 1e-5) {
ci <- ci + 1; if (ci > 5) break
FF <- fy(ytm, cf, tm)
dFF <- (fy(ytm + h, cf, tm) - FF) / h
dF <- FF / dFF
ytm <- ytm - dF
}
if (ytm < 0)
ytm <- 0.99
ytm
}
###################################################
### code chunk number 217: NMOFman.Rnw:4612-4638
###################################################
OF3 <- function(param, Data) {
tm <- Data$tm
rM <- Data$rM
cfMatrix<- Data$cfMatrix
nB <- dim(cfMatrix)[2L]
zrates <- Data$model(param,tm); aux <- 1e10
if ( all(zrates > 0,
!is.na(zrates))
) {
diFa <- 1 / ((1 + zrates/100)^tm)
b <- diFa %*% cfMatrix
r <- numeric(nB)
if ( all(!is.na(b),
diFa < 1,
diFa > 0,
b > 1)
) {
for (bb in 1:nB) {
r[bb] <- compYield(c(-b[bb], cfMatrix[ ,bb]), c(0,tm))
}
aux <- abs(r - rM)
aux <- sum(aux)
}
}
aux
}
###################################################
### code chunk number 218: NMOFman.Rnw:4648-4653
###################################################
betaTRUE <- c(5,-2,1,10,1,3)
yM <- NSS(betaTRUE, tm)
diFa <- 1 / ( (1 + yM/100)^tm )
bM <- diFa %*% cfMatrix
rM <- apply(rbind(-bM, cfMatrix), 2, compYield, c(0, tm))
###################################################
### code chunk number 219: NMOFman.Rnw:4658-4678
###################################################
Data <- list(rM = rM, tm = tm,
cfMatrix = cfMatrix,
model = NSS,
min = c( 0,-15,-30,-30,0 ,2.5),
max = c(15, 30, 30, 30,2.5,5 ),
ww = 0.1,
fy = fy)
algo <- list(nP = 100L,
nG = 1000L,
F = 0.50,
CR = 0.99,
min = c( 0,-15,-30,-30,0 ,2.5),
max = c(15, 30, 30, 30,2.5,5 ),
pen = penalty,
repair = NULL,
loopOF = TRUE,
loopPen = FALSE,
loopRepair = FALSE,
printBar = FALSE,
printDetail = FALSE)
###################################################
### code chunk number 220: NMOFman.Rnw:4681-4684
###################################################
sol <- DEopt(OF = OF3, algo = algo, Data = Data)
max(abs(Data$model(sol$xbest,tm) - Data$model(betaTRUE,tm)))
sol$OFvalue
###################################################
### code chunk number 221: NMOFman.Rnw:4689-4697
###################################################
s0 <- algo$min + (algo$max - algo$min) * runif(length(algo$min))
sol2 <- nlminb(s0, OF3, Data = Data,
lower = algo$min,
upper = algo$max,
control = list(eval.max = 50000L,
iter.max = 50000L))
max(abs(Data$model(sol2$par,tm) - Data$model(betaTRUE,tm)))
sol2$objective
###################################################
### code chunk number 222: NMOFman.Rnw:4700-4709
###################################################
par(ps = 8, bty = "n", las = 1, tck = 0.01,
mgp = c(3, 0.5, 0), mar = c(4, 4, 1, 1))
plot(tm, yM, xlab = "maturities in years", ylab = "yields in %")
lines(tm,Data$model(sol$xbest,tm), col = "blue")
lines(tm,Data$model(sol2$par,tm), col = "darkgreen", lty = 2)
legend(x = "bottom", legend = c("true yields","DE","nlminb"),
col = c("black", "blue", "darkgreen"),
pch = c(1, NA, NA), lty = c(0,1,2))
###################################################
### code chunk number 223: NMOFman.Rnw:4713-4715
###################################################
betaTRUE
round(sol$xbest,3)
###################################################
### code chunk number 224: NMOFman.Rnw:4728-4729
###################################################
rm(list = ls())
###################################################
### code chunk number 225: NMOFman.Rnw:4733-4735
###################################################
require("NMOF")
set.seed(94679)
###################################################
### code chunk number 226: NMOFman.Rnw:4757-4776
###################################################
randomData <- function(p = 200L, ## number of available regressors
n = 200L, ## number of observations
maxReg = 10L, ## max. number of included regressors
s = 1, ## standard deviation of residuals
constant = TRUE ) {
X <- array(rnorm(n * p), dim = c(n, p))
if (constant)
X[ ,1L] <- 1
k <- sample.int(maxReg, 1L) ## the number of true regressors
K <- sort(sample.int(p, k)) ## the set of true regressors
betatrue <- rnorm(k) ## the true coefficients
## the response variable y
y <- X[ ,K] %*% as.matrix(betatrue) + rnorm(n, sd = s)
list(X = X, y = y, betatrue = betatrue, K = K, n = n, p = p)
}
###################################################
### code chunk number 227: NMOFman.Rnw:4780-4782
###################################################
rD <- randomData(p = 100L, n = 200L, s = 1,
constant = TRUE, maxReg = 10L)
###################################################
### code chunk number 228: NMOFman.Rnw:4785-4791
###################################################
Data <- list(X = rD$X,
y = rD$y,
n = rD$n,
p = rD$p,
maxk = 30L, ## maximum number of regressors included in model
lognn = log(rD$n)/rD$n)
###################################################
### code chunk number 229: NMOFman.Rnw:4794-4798
###################################################
x0 <- logical(Data$p)
temp <- sample.int(Data$maxk, 1L)
temp <- sample.int(Data$p, temp)
x0[temp] <- TRUE
###################################################
### code chunk number 230: NMOFman.Rnw:4807-4808
###################################################
rD$K
###################################################
### code chunk number 231: NMOFman.Rnw:4811-4812
###################################################
which(x0)
###################################################
### code chunk number 232: NMOFman.Rnw:4825-4829
###################################################
result1 <- lm(Data$y ~ -1 + Data$X[ ,x0])
result2 <- qr.solve(Data$X[ ,x0], Data$y)
## ... coefficients should be the same
all.equal(as.numeric(coef(result1)), as.numeric(result2))
###################################################
### code chunk number 233: NMOFman.Rnw:4832-4838
###################################################
require("rbenchmark")
benchmark(lm(Data$y ~ -1 + Data$X[ ,x0]),
qr.solve(Data$X[ ,x0], Data$y),
columns = c("test", "elapsed", "relative"),
order = "test",
replications = 1000L)
###################################################
### code chunk number 234: NMOFman.Rnw:4851-4856
###################################################
OF <- function(x, Data) {
q <- qr(Data$X[ ,x])
e <- qr.resid(q, Data$y)
log(crossprod(e)/Data$n) + sum(x) * Data$lognn
}
###################################################
### code chunk number 235: NMOFman.Rnw:4859-4860
###################################################
OF(x0, Data)
###################################################
### code chunk number 236: NMOFman.Rnw:4866-4874
###################################################
neighbour <- function(xc, Data) {
xn <- xc
ex <- sample.int(Data$p, 1L)
xn[ex] <- !xn[ex]
sumx <- sum(xn)
if ( sumx < 1L || (sumx > Data$maxk) )
xc else xn
}
###################################################
### code chunk number 237: NMOFman.Rnw:4878-4881
###################################################
OF(neighbour(x0, Data), Data)
OF(neighbour(x0, Data), Data)
OF(neighbour(x0, Data), Data)
###################################################
### code chunk number 238: NMOFman.Rnw:4886-4893
###################################################
algo <- list(nT = 10L, ## number of thresholds
nS = 200L, ## number of steps per threshold
nD = 1000L, ## number of random steps to compute thresholds
neighbour = neighbour,
x0 = x0,
printBar = FALSE)
system.time(sol1 <- TAopt(OF, algo = algo, Data = Data))
###################################################
### code chunk number 239: NMOFman.Rnw:4899-4902
###################################################
sol1$OFvalue
which(sol1$xbest) ## the selected regressors
rD$K ## the true regressors
###################################################
### code chunk number 240: NMOFman.Rnw:4909-4913
###################################################
xtrue <- logical(Data$p)
xtrue[rD$K] <- TRUE
OF(sol1$xbest, Data)
OF(xtrue, Data)
###################################################
### code chunk number 241: NMOFman.Rnw:4920-4929
###################################################
restarts <- 50L
algo$printDetail <- FALSE
res <- restartOpt(TAopt, n = restarts,
OF = OF, algo = algo, Data = Data,
method = "snow", cl = 2)
par(bty = "n", las = 1,mar = c(3,4,0,0),
ps = 8, tck = 0.001, mgp = c(3, 0.5, 0))
plot(ecdf(sapply(res, `[[`, "OFvalue")), ## extract solution quality
cex = 0.4, main = "", ylab = "", xlab = "")
###################################################
### code chunk number 242: NMOFman.Rnw:4933-4941
###################################################
xbestAll <- sapply(res, `[[`, "xbest") ## extract all solutions
inclReg <- which(rowSums(xbestAll) > 0L) ## get included regressors
inclReg <- sort(union(rD$K, inclReg))
data.frame(regressor = inclReg,
`included` = paste(rowSums(xbestAll)[inclReg], "/",
restarts, sep = ""),
`true regressor?` = inclReg %in% rD$K,
check.names = FALSE)
###################################################
### code chunk number 243: NMOFman.Rnw:4963-5000
###################################################
dim(rD$X)
neighbour2 <- function(xc, Data) {
if ((sumx <- sum(x0)) >= Data$maxk)
ex <- sample(which(x0), 1L)
else if (sumx == 1L)
ex <- sample(which(!x0), 1L)
else
ex <- sample.int(Data$p, 1L)
xc[ex] <- !xc[ex]
xc
}
neighbour <- function(xc, Data) {
xn <- xc
ex <- sample.int(Data$p, 1L)
xn[ex] <- !xn[ex]
sumx <- sum(xn)
if ( sumx < 1L || (sumx > Data$maxk) )
xc else xn
}
algo <- list(nT = 20L, ## number of thresholds
nS = 200L, ## number of steps per threshold
nD = 1000L, ## number of random steps to compute thresholds
neighbour = neighbour,
x0 = x0, q= 0.5,
printBar = FALSE)
system.time(sol1 <- TAopt(OF, algo = algo, Data = Data))
plot(cummin(sol1$Fmat[ ,2L]), type = "l", log = "y")
##rD
OF <- function(x, Data) {
q <- qr(Data$X[ ,x])
e <- qr.resid(q, Data$y)
crossprod(e)
}
###################################################
### code chunk number 244: NMOFman.Rnw:5022-5029
###################################################
bsm <- function(S, X, tau, r, q, vol, I = 1) {
d1 <- (log(S/X) + (r - q + vol^2/2) * tau)/(vol * sqrt(tau))
d2 <- d1 - vol * sqrt(tau)
list(value = I * (S * exp(-q * tau) * pnorm(I * d1) -
X * exp(-r * tau) * pnorm(I * d2)),
vega = S * exp(-q*tau) * dnorm(d1 * I) * sqrt(tau))
}
###################################################
### code chunk number 245: NMOFman.Rnw:5035-5043
###################################################
S <- 99 ## sport
X <- 100 ## strike
r <- 0.01
q <- 0.0
tau <- 0.25
vol <- 0.2
I <- 1 ## a call (-1 for a put)
unlist(bsm(S, X, tau, r, q, vol, I))
###################################################
### code chunk number 246: NMOFman.Rnw:5051-5056
###################################################
tmp <- unlist(vanillaOptionEuropean(S = S, X = X, tau = tau,
r = r, q = q, v = vol^2,
type =
ifelse(I == 1, "call", "put")))
tmp[c("value","vega")]
###################################################
### code chunk number 247: NMOFman.Rnw:5064-5072
###################################################
S <- 99
X <- 100
r <- 0.01
q <- 0.01
tau <- 0.1
I <- 1
vol <- 0.247
(price <- bsm(S, X, tau, r, q, vol, I)$value)
###################################################
### code chunk number 248: NMOFman.Rnw:5080-5092
###################################################
impliedVol <- function(price, S, X, tau, r, q, vol0 = 0.15, I = 1,
tol = 1e-4, maxit = 10) {
for (i in seq_len(maxit)) {
tmp <- bsm(S, X, tau, r, q, vol0, I)
step <- (tmp$value - price)/tmp$vega
vol0 <- vol0 - step
if (all(abs(step) < tol))
break
}
vol0
}
###################################################
### code chunk number 249: NMOFman.Rnw:5099-5100
###################################################
impliedVol(price, S, X, tau, r, q, vol, I)
###################################################
### code chunk number 250: NMOFman.Rnw:5105-5107
###################################################
vanillaOptionImpliedVol(exercise = "european", price, S, X, tau, r,
q, type = "call")
###################################################
### code chunk number 251: NMOFman.Rnw:5114-5121 (eval = FALSE)
###################################################
## benchmark(iV = impliedVol(price, S, X, tau, r, q, runif(1L) + 0.01,I),
## vanOptIV = vanillaOptionImpliedVol(exercise = "european",
## price, S, X, tau, r,
## q, tauD = 0, D = 0, type = "call",
## M = 101, uniroot.info = FALSE),
## columns = c("test", "elapsed", "relative"),
## replications = 1e2, order = "relative")
###################################################
### code chunk number 252: NMOFman.Rnw:5131-5136
###################################################
S <- rep(99, 21) ## spot
X <- 90:110 ## strike
r <- 0.01; q <- 0.02
tau <- 0.2; vol <- 0.24; I <- 1
data.frame(S = S, X = X, bsm(S, X, tau, r, q, vol, I))
###################################################
### code chunk number 253: NMOFman.Rnw:5142-5145
###################################################
Xvec <- 80:120
tauvec <- c(c(3, 6, 9)/12, ## 3, 6, 9 months
1, 2, 3, 4, 5) ## 1..5 years
###################################################
### code chunk number 254: NMOFman.Rnw:5151-5158
###################################################
loop <- function() {
callprices <- array(NA, dim = c(length(Xvec), length(tauvec)))
for (X in Xvec)
for (tau in tauvec)
callprices[X == Xvec, tau == tauvec] <- bsm(S,X,tau,r,q,vol)$value
callprices
}
###################################################
### code chunk number 255: NMOFman.Rnw:5163-5169
###################################################
vect <- function() {
tmp <- expand.grid(Xvec,tauvec)
callprices <- bsm(S, tmp[[1L]], tmp[[2L]], r, q, vol, I)$value
dim(callprices) <- c(length(Xvec), length(tauvec))
callprices
}
###################################################
### code chunk number 256: NMOFman.Rnw:5174-5185
###################################################
S <- 101
Xvec <- 80:120
tauvec <- c(c(3, 6, 9)/12, ## 3, 6, 9 months
1, 2, 3, 4, 5) ## 1..5 years
r <- 0.01; q <- 0.01
tau <- 0.25; vol <- 0.2; I <- 1
callprices1 <- loop()
callprices2 <- vect()
all.equal(callprices1, callprices2)
###################################################
### code chunk number 257: NMOFman.Rnw:5191-5194 (eval = FALSE)
###################################################
## benchmark(loop(), vect(),
## columns = c("test", "elapsed", "relative"),
## replications = 1e2, order = "relative")
###################################################
### code chunk number 258: NMOFman.Rnw:5201-5212
###################################################
S <- rep(99,21) ## spot
X <- 90:110
r <- 0.01
q <- 0.02
tau <- runif(21)
vol <- (runif(21)+0.2)/3
ivol <- impliedVol(bsm(S, X, tau, r, q, vol, I)$value,
S, X, tau, r, q, vol = 0.2,
I, tol = 1e-5, maxit = 5)
data.frame(S = S, X = X, vol = vol, ivol = ivol, diff = abs(vol-ivol))
###################################################
### code chunk number 259: NMOFman.Rnw:5219-5223 (eval = FALSE)
###################################################
## benchmark(iv = impliedVol(bsm(S, X, tau, r, q, vol, I)$value,
## S, X, tau, r, q, tol = 1e-5, maxit = 5),
## columns = c("test", "elapsed"),
## replications = 1e3, order = "relative")
###################################################
### code chunk number 260: NMOFman.Rnw:5319-5345 (eval = FALSE)
###################################################
## valDate <- as.Date("2012-02-10")
##
## calls <- optionData$pricesCall
## puts <- optionData$pricesPut
## spot <- optionData$index
##
## ## strikes
## strikes <- as.numeric(row.names(calls))
##
## ## compute interest rates
## expDates <- as.Date(paste0(colnames(calls),"15"),
## format = "%Y%m%d")
## tau <- as.numeric((expDates - valDate)/365)
## r <- NSS(optionData$NSSpar, tau)
##
## ## compute basis
## notNA <- !is.na(calls) & !is.na(puts)
## i <- 2
## q <- numeric(ncol(calls))
## for (i in seq_along(q))
## q[i] <- mean(-log((calls[,i] - puts[,i] +
## exp(-r[i]*tau[i]) * strikes)/spot)/tau[i],
## na.rm = TRUE)
##
## ## check with future
## ((log(spot)-log(optionData$future[2]))+r[2]*tau[2])/tau[2]
###################################################
### code chunk number 261: NMOFman.Rnw:5349-5363 (eval = FALSE)
###################################################
## pricevec <- as.numeric(calls)
## notNA <- !is.na(pricevec)
## Xvec <- rep(strikes, ncol(calls))
## tauvec <- rep(tau, each = nrow(calls))
## rvec <- rep(r, each = nrow(calls))
## qvec <- rep(q, each = nrow(calls))
## ivol <- impliedVol(pricevec[notNA],
## spot, Xvec[notNA], tauvec[notNA],
## rvec[notNA], qvec[notNA], vol = 0.25, maxit = 5)
##
## vols <- numeric(length(pricevec))
## vols[notNA] <- ivol
## data.frame(pricevec[notNA],
## spot, Xvec[notNA], tauvec[notNA], vols=ivol)
###################################################
### code chunk number 262: NMOFman.Rnw:5372-5391
###################################################
S <- 100 ## spot
X <- 100 ## strike
tau <- 1 ## time-to-maturity
r <- 0.02 ## interest rate
q <- 0.02 ## dividend rate
v <- 0.2 ## volatility
## the closed-form solution
callBSM <- function(S,X,tau,r,q,v) {
d1 <- (log(S/X) + (r - q + v^2 / 2)*tau) / (v*sqrt(tau))
d2 <- d1 - v*sqrt(tau)
S * exp(-q * tau) * pnorm(d1) - X * exp(-r * tau) * pnorm(d2)
}
callBSM(S,X,tau,r,q,v)
## with the characteristic function
callCF(cf = cfBSM, S = S, X = X, tau = tau, r = r, q = q,
v = v^2, ## variance, not vol
implVol = TRUE)
###################################################
### code chunk number 263: NMOFman.Rnw:5397-5400
###################################################
Xvec <- 80:120
tauvec <- c(c(3, 6, 9)/12, ## 3, 6, 9 months
1, 2, 3, 4, 5) ## 1..5 years
###################################################
### code chunk number 264: NMOFman.Rnw:5404-5411
###################################################
fun1 <- function() {
callprices <- array(NA, dim = c(length(Xvec), length(tauvec)))
for (X in Xvec)
for (tau in tauvec)
callprices[X == Xvec, tau == tauvec] <- callBSM(S,X,tau,r,q,v)
callprices
}
###################################################
### code chunk number 265: NMOFman.Rnw:5415-5421
###################################################
fun2 <- function() {
tmp <- expand.grid(Xvec,tauvec)
callprices <- callBSM(S, tmp[[1]], tmp[[2]], r, q, v)
dim(callprices) <- c(length(Xvec), length(tauvec))
callprices
}
###################################################
### code chunk number 266: NMOFman.Rnw:5424-5428
###################################################
callprices1 <- fun1()
callprices2 <- fun2()
all.equal(callprices1, callprices2)
###################################################
### code chunk number 267: NMOFman.Rnw:5432-5438
###################################################
system.time(
for (i in 1:100)
ignore <- fun1() )
system.time(
for (i in 1:100)
ignore <- fun2() )
###################################################
### code chunk number 268: NMOFman.Rnw:5443-5491
###################################################
priceMatrix <- function(cf, S, Xvec, tauvec, r, q = 0, ...,
nodes = NULL, weights = NULL, n = 200) {
if (is.null(nodes)) {
tmp <- xwGauss(n)
tmp <- changeInterval(tmp$nodes, tmp$weights,
oldmin = -1, oldmax = 1,
newmin = 0, newmax = 200)
nodes <- tmp$nodes
weights <- tmp$weights
}
callprices <- array(NA, dim = c(length(Xvec), length(tauvec)))
tmpmat <- array(NA, dim = c(length(Xvec), length(weights)))
inodes <- 1i * nodes
itau <- 0L
for (tau in tauvec) {
itau <- itau + 1L
cfi <- S * exp((r - q) * tau)
cf1 <- cf(nodes - 1i, S, tau, r, q, ...)/inodes/cfi
cf2 <- cf(nodes, S, tau, r, q, ...)/inodes
iX <- 0L
for (X in Xvec) {
iX <- iX + 1L
if (itau == 1L)
tmpmat[iX, ] <- exp(-inodes * log(X))
P1 <- 0.5 + weights %*% Re(tmpmat[iX, ] * cf1)/pi
P2 <- 0.5 + weights %*% Re(tmpmat[iX, ] * cf2)/pi
callprices[iX, itau] <-
exp(-q*tau) * S * P1 - exp(-r*tau) * X * P2
}
}
callprices
}
v <- .2
system.time(
for (i in 1:100)
ignore <- priceMatrix(cf = cfBSM, S, Xvec, tauvec, r, q = q,
v = 0.2^2, n = 50) )
require("compiler")
priceMatrix2 <- cmpfun(priceMatrix)
system.time(
for (i in 1:100)
ignore <- priceMatrix2(cf = cfBSM, S, Xvec, tauvec, r, q = q,
v = 0.2^2, n = 50) )
###################################################
### code chunk number 269: NMOFman.Rnw:5494-5524
###################################################
cfp <- priceMatrix(cf = cfBSM, S, Xvec, tauvec, r, q = q,
v = 0.2^2, n = 100)
callprices1[1:5, 1:5]
callprices2[1:5, 1:5]
cfp[1:5, 1:5]
all.equal(callprices1, cfp)
system.time(
for (i in 1:100)
ignore <- fun1() )
system.time(
for (i in 1:100)
ignore <- fun2() )
system.time(
for (i in 1:100)
ignore <- priceMatrix(cf = cfBSM, S, Xvec, tauvec, r, q = q,
v = 0.2^2, n = 50) )
tmp <- xwGauss(50)
tmp <- changeInterval(tmp$nodes, tmp$weights,
oldmin = -1, oldmax = 1,
newmin = 0, newmax = 200)
nodes <- tmp$nodes
weights <- tmp$weights
system.time(
for (i in 1:100)
ignore <- priceMatrix(cf = cfBSM, S, Xvec, tauvec, r, q = q,
v = 0.2^2,
nodes = nodes, weights = weights) )
###################################################
### code chunk number 270: NMOFman.Rnw:5612-5614 (eval = FALSE)
###################################################
## install.packages("NMOF") ## CRAN
## install.packages("NMOF", repos = "http://enricoschumann.net/R")
###################################################
### code chunk number 271: NMOFman.Rnw:5622-5624 (eval = FALSE)
###################################################
## require("NMOF")
## showExample("exampleOF.R")
###################################################
### code chunk number 272: NMOFman.Rnw:5655-5658 (eval = FALSE)
###################################################
## require("utils")
## bug.report("[NMOF] Unexpected behaviour in function XXX",
## maintainer("NMOF"), package = "NMOF")
###################################################
### code chunk number 273: NMOFman.Rnw:5666-5667
###################################################
toLatex(sessionInfo())
Cannot open file '': No such file or directory
Warning message:
In file(con, "r") :
file("") only supports open = "w+" and open = "w+b": using the former
Number of unit tests: 0
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.