Nothing
# Code implementing different log optimal portfolio selection
# algorithms, especially linked to universal portfolios
#
require(xts)
require(quadprog)
require(FNN)
# check for compatible dimensions
# b is either a vector or a matrix of weights
# x is n (timesteps) * m (assets)
# we need b = 1 * m or m * 1 or n * m
is.compatible <- function (b, x) {
if(is.vector(b)) { return(dim(x)[2] == length(b)) }
if(is.array(b)) { return(all(dim(b) == dim(x))) }
return(FALSE)
}
# sum of weights must be 1 per line
normalize <- function (b) {
if (is.vector(b)) { return (b/sum(b)) }
if (is.array(b)) { return (b/apply(b,1,sum)) }
return(NULL)
}
# take a wealth/price sequence and calculate the wealth/price relative
w2x <- function (w) {
if (is.xts(w)) {
x <- w / lag(w)
x[1,] <- w[1,]
return(x)
}
x <- w
if (is.vector(w)) {
n <- length(w)
x[2:n] <- w[2:n] / w[1:(n-1)]
return(x)
}
if (is.array(w)) {
n <- dim(w)[1]
x[2:n,] <- w[2:n,] / w[1:(n-1),]
return(x)
}
return(NULL)
}
# take a price/wealth relative relative and calculate the price/wealth
x2w <- function (x) {
if (is.xts(x)) { return(cumprod(x)) }
if (is.array(x)) { return(apply(x,2,cumprod)) }
if (is.vector(x)) { return(cumprod(x)) }
return(NULL)
}
# crp and bh are very close, so we implement them as the same procedure
# bh works on w, crp on x, but otherwise they are similar
# - b can be either a vector (constant) or a matrix (weight at each instant in time)
# if b is NULL, we use an uniform vector loading, simple reference
# - x is a matrix of price relative, normally an xts
crp_bh <- function (x, b = NULL, bh = FALSE) {
if (is.null(b)) { b <- rep(1/dim(x)[2],dim(x)[2]) }
if(!is.compatible(b,x)) { stop("b and x not compatible in crp(x,b)", Call.=TRUE) }
b <- normalize(b)
if (bh) { x <- x2w(x) }
if (is.vector(b)) { xp <- x %*% b }
else { xp <- apply(b*x, 1, sum) }
if (is.xts(x)) { xp <- xts(xp, order.by=index(x)) }
if (bh) { return(xp) }
else { return(x2w(xp)) }
}
crp <- function (x, b = NULL) {
return (crp_bh(x, b, bh = FALSE))
}
bh <- function (x, b = NULL) {
return (crp_bh(x, b, bh = TRUE))
}
oracle <- function (x) { # non causal maximum gain
if(is.vector(x)) { return(x2w(x)) }
xp <- apply(x,1,max)
if (is.xts(x)) { return(xts(x2w(xp), order.by=index(x))) }
return(x2w(xp))
}
best.asset <- function (x) {
w <- x2w(x)
if(is.vector(x)) { return(w) }
iBest <- which.max(w[dim(x)[1],])
return(w[,iBest])
}
worst.asset <- function (x) {
w <- x2w(x)
if(is.vector(x)) { return(w) }
iWorst <- which.min(w[dim(x)[1],])
return(w[,iWorst])
}
best.envelope <- function (x) {
w <- x2w(x)
if(is.vector(x)) { return(w) }
wm <- apply(w,1,max)
if(is.xts(x)) { return(xts(wm, order.by=index(x))) }
return(wm)
}
worst.envelope <- function (x) {
w <- x2w(x)
if(is.vector(x)) { return(w) }
wm <- apply(w,1,min)
if(is.xts(x)) { return(xts(wm, order.by=index(x))) }
return(wm)
}
round.to.zero <- function (b, do.round, eps) {
if (do.round) { b[abs(b) < eps] <- 0 }
return(normalize(b))
}
# use quadratic programming approximation for fast results in finding best crp
# this is sometimes called semilog
# return the optimal vector of weights
bcrp.qp <- function (x, clean.up = TRUE, clean.up.eps = 1e-10 ) {
if (is.vector(x)) { return(c(1)) }
Dmat <- t(x-1) %*% (x-1)
dvec <- apply(x-1,2,sum)
lb <- dim(x)[[2]]
Amat <- matrix(0,lb,lb+1)
Amat[,1] <- 1
for (i in 1:lb) { Amat[i,i+1] <- 1 }
bvec <- 0 * 1:(lb+1)
bvec[1] <- 1
b <- solve.QP(Dmat, dvec, Amat, bvec, meq=1)
bc <- round.to.zero(b$solution, clean.up, clean.up.eps)
if(any(is.na(bc))) { return(b$solution) }
else { return(bc) }
}
# transform axis to insure non negative values, work except that it is not possible to reach 0 values
# alternate approach would use constrained optimization
CRPoptrescale <- function (b, x) {
b <- exp(b)
return(-crp(x,b)[dim(x)[1]])
}
# find the optimal CRP weights using optim (can be slow)
# optim is unconstrained, and we use a trick to map any value as positive
# then to remap it into the simplex, see CRPoptrescale
bcrp.optim <- function (x, maxit=20, clean.up = TRUE, clean.up.eps = 1e-10, fast.only = FALSE ) {
if(is.vector(x)) { return(c(1)) }
if (dim(x)[1] > dim(x)[2]) { b0 <- try(bcrp.qp(x)) }
else { b0 <- rep(1,dim(x)[2]); fast.only = FALSE }
if (inherits(b0,"try-error")) { # bcrp.qp may fail for ill conditioned problems, use best stock as start
b0 <- rep(0,dim(x)[2])
b0[which.max(apply(x,2,prod))] <- 1
fast.only = FALSE
}
if (!fast.only) {
solution <- optim(log(b0+clean.up.eps), CRPoptrescale, gr=NULL, method="BFGS", lower=-Inf, upper=Inf,control = list(maxit=maxit), hessian=FALSE, x)
b0 <- exp(solution$par)
}
return(round.to.zero(normalize(b0), clean.up, clean.up.eps))
}
# use a "urns and balls" approach, a is a vector of balls per urn
# return NULL when done (all balls in the rightmost urn)
next.simplex.point <- function (a) {
nurns <- length(a)
for (i in 1:nurns) {
if (a[i] > 0) {
if (i == nurns) { return(NULL) }
else {
a[i+1] <- 1 + a[i+1]
a[1] <- a[i] - 1
if (i != 1) { a[i] <- 0 }
return(a)
}
}
}
}
# n urns and m balls, i.e. number of simplex points that will be sampled
count.grid.points <- function (n, m) {
pts <- 1
if (n < 2) { return(1) }
for (i in 1:(n-1)) { pts <- pts * (m+i) / i }
return(pts)
}
# n is the number of sample intervals per dimension, i.e. step is 1/n
universal.cover <- function (x, n) {
if (is.vector(x)) { return(x2w(x)) }
m <- dim(x)[2]
b <- 0 * 1:m
b[1] <- n
npts <- 0
w <- 0 * x[,1]
repeat {
w <- w + crp(x, b)
npts <- npts + 1
b <- next.simplex.point(b)
if (is.null(b)) {
return(w/npts)
}
}
}
# Ishijima methods for uniform and Dirichlet sampling of the simplex
# uniform used the direct method (with sort), not the sequential method with exponent (could be done in fact)
# n is the number of samples
# m is the number of dimensions
uniform.simplex.sampling <- function (n,m) {
if(m==1) { return(runif(n)) }
x <- array(runif((m-1) * n), dim=c(n,(m-1))) # prepare for the operation of apply
if (m > 2) {
x <- t(apply(x,1,sort))
y <- apply(x,1,diff)
if (m>3) { y <- t(y) }
b <- cbind(x[,1], y, x[,m-1])
return(b)
}
return (cbind(x, 1-x))
}
# n is the number of samples
# m is the number of dimensions
# alpha is the parameter for the gamma
dirichlet.simplex.sampling <- function (n, m, alpha=0.5) {
if(m==1) { return (rgamma(n, alpha)) }
x <- array(rgamma(m * n, alpha), dim=c(n,m))
sm <- x %*% rep(1, m)
x/as.vector(sm)
}
# method is "uniform" or "dirichlet"
# generation of points in simplex with correct density is based on Ishijima article
# implementation for dirichlet is however taken from rdirichlet
# n is the number of sample points
universal.cover.random <- function (x, n, method) { # should also return a vector of combined b
if(is.vector(x)) { return(w2x(x)) }
m <- dim(x)[2]
w <- 0 * x[,1]
if (method == "uniform") { bs <- uniform.simplex.sampling(n,m) }
if (method == "dirichlet") { bs <- dirichlet.simplex.sampling(n,m) }
for (i in 1:n) { w <- w + crp(x, bs[i,]) }
return (w/n)
}
# multiplicative updates algorithm, does not match published results
# supported methods are: "EG" and "chi2"
# "exact" is planned
# very simple algorithm, but require a loop
mult.upgrade <- function (x, eta, method) {
if(is.vector(x)) { return(c(1)) }
if(is.xts(x)) { x <- array(x,dim=dim(x)) }
b <- 0 * x
b[1,] <- 1/dim(x)[2]
if (method == "EG") {
for (t in 2:dim(x)[1]) {
b[t,] <- b[t-1,] * exp(eta * x[t-1,] / (b[t-1,] %*% x[t-1,]))
b[t,] <- b[t,] / sum(b[t,])
}
return(b)
}
if (method == "chi2") {
for (t in 2:dim(x)[1]) {
b[t,] <- b[t-1,] * (eta * (x[t-1,] / (b[t-1,] %*% x[t-1,]) - 1) + 1)
b[t,] <- b[t,] / sum(b[t,])
}
return(b)
}
if (method == "exact") {
print("Method 'exact' is not currently supported in function mult.upgrade")
return(NULL)
}
return(NULL)
}
# Yoram Singer switching portfolio
# method is "fixed" or "adaptive", "adaptive" does not reproduce published results
switching.portfolio <- function(x, gamma, method) {
# w[t,i] store the value of asset i after application of x[t,i]
w <- array(0, dim=dim(x))
nAssets <- dim(x)[2]
w[1,] <- 1/nAssets * x[1,]
if (method == "fixed") {
for (t in 2:dim(x)[1]) {
w[t,] <- ((1 - gamma - gamma / (nAssets-1)) *
w[t-1,] +
(gamma / (nAssets-1)) * sum(w[t-1,])) * x[t,]
}
w <- apply(w,1,sum)
if(is.xts(x)) { return(xts(w, order.by=index(x))) }
return(w)
}
if (method == "adaptive") {
# wtt0[t0,i] store the value of the pure strategy starting at t0 before redistribution
# ghdt stands for gamma hat delta t
wtt0 <- w
for (t in 2:dim(x)[1]) {
past.index <- 1:(t-1)
ghdt <- 0.5 / (1 + seq(t-1, 1, by = -1))
out <- wtt0[past.index,] * ghdt # rely on recycling
if (t == 2) { out.by.asset <- out/(nAssets-1) }
else { out.by.asset <- apply(out,2,sum) / (nAssets-1) }
wtt0[t,] <- -out.by.asset + sum(out.by.asset)
wtt0[past.index,] <- wtt0[past.index,] - out
wtt0[1:t,] <- wtt0[1:t,] * rep(x[t,], each=t)
w[t,] <- apply(wtt0[1:t,],2,sum)
}
w <- apply(w,1,sum)
if(is.xts(x)) { return(xts(w, order.by=index(x))) }
return(w)
}
return(NULL)
}
# roll.bcrp returns the BCRP calculated at a number of points
# and taking into account a window (potentially infinite)
# wl is the window length
# return is a list with
# - the indices of the calculated starts
# - the indices of the calculated ends
# - a matrix of weights
roll.bcrp <- function(x, from=1, by=1, wl=NULL, fast.only=TRUE) {
# smallish bug, need at least 2 rows for correct operation
n <- dim(x)[1]
m <- dim(x)[2]
from <- max(2,from)
ends <- seq(from,n,by)
nw <- length(ends)
b <- array(1/m,dim=c(nw,m))
if(is.null(wl)) { wl <- n }
starts <- pmax(ends-wl+1,1)
for (i in 1:nw) {
b[i,] <- bcrp.optim(x[starts[i]:ends[i],], fast.only=fast.only)
}
return(list(starts, ends, b))
}
# EMA on the Successive CRP to avoid some instability
wscrp <- function (x, from=1, by=1, alpha=0.99, fast.only=TRUE) {
m <- dim(x)[2]
n <- dim(x)[1]
b <- 0 * x + 1/m
from <- max(from,2)
i <- from+1
# smallish bug, need at least 2 rows for bcrp.optim
# the problem is that bcrp.optim has been rewritten to assume that a vector is
# a one column vector, while here it would be a one row vector
# unfortunately a one row slice becomes a vector and cannot be dsitinguished from
# a one column vector. Solution would be to force the slice to be an array of the
# correct dimension, but this is annoying. Note that this is only the case for
# arrays, xts does distinguish between the two cases.
b0 <- b[1,]
while(i < n) {
b0 <- bcrp.optim(x[1:(i-1),], fast.only=fast.only) * (1-alpha) + b0 * alpha
b0[is.nan(b0)] <- 1/m
j <- min(n,i+by)
bt <- array(rep(b0,j-i+1),dim=c(m,j-i+1))
b[i:j,] <- t(bt)
i <- j
}
return(crp(x, b))
}
scrp <- function (x, from=1, by=1, fast.only=TRUE) {
return(wscrp(x, from=from, by=by, alpha=0, fast.only=fast.only))
}
# nearest neighbor approach, from Gyorfi, only calculate the porftolio
# themselves at this time, not the weights
# x contains the price relative
# pls is either
# - a scalar, in which case it is the value L in Gyorfi and the pl sequence is 0.02 + 0.05 * (l-1)/(L-1), l in 1:L
# - directly the pl sequence
# ks is either
# - a scalar, in which case it is the value K in Gyorfi and the k sequence is 1:K
# - directly the k sequence
# qkl is the a priori loading on the K*L set of experts, if NULL uniform loading is used
#
nn.gyorfi <- function(x, pls, ks, qkl=NULL) {
if (is.xts(x)) { x <- array(x,dim=dim(x)) }
if (is.array(pls)) { L <- length(pls) }
else { L <- pls ; pls <- 0.02 + 0.5 / (L-1) * ((1:L)-1) }
if (is.array(ks)) { K <- length(ks) }
else { K <- ks ; ks <- 1:K }
if(is.null(qkl)) { qkl <- array(1/K/L, dim=c(K,L)) }
else { # check for compatibility
if(any(dim(qkl) != c(K,L))) {
print("qkl loading is not compatible with array of experts")
return(NULL)
}
qkl <- qkl/sum(qkl)
}
nSamples <- dim(x)[1]
nAssets <- dim(x)[2]
b <- array(1/nAssets,dim=c(K,L,nSamples,nAssets))
w <- array(0,dim=c(K,L,nSamples))
# the order of the loop is
# - k first, to construct the target array
# - then i (time), to calculate the set of nearest neighbors
# - then l to extract a subset of the nearest neighbors and optimize
for (ik in seq_along(ks)) {
k <- ks[ik]
# construct the target sequence by column juxtaposition
for (lag in 0:(k-1)) {
if (lag == 0) { xnn <- x[1:(nSamples-k+1),] }
else { xnn <- cbind(xnn, x[(1+lag):(nSamples-k+1+lag),]) }
}
# initial part uses uniform, could extend to more than k
for (t in (k+2):(nSamples-k-1)) {
# calculate the set of number of nearest neighbors and their maximum
# extract the maximum number of nearest neighbor
ls <- floor((t-k+1) * pls)
if (max(ls) < 1) { next }
nnsk <- knnx.index(xnn[1:(t-k),],array(xnn[t-k+1,],dim=c(1,nAssets*k)),max(ls),algorithm="kd_tree")
for (il in seq_along(ls)) {
l <- ls[il]
if (l < 1) { next }
xcrp <- x[nnsk[1:l]+k,]
b[ik,il,t+1,] <- bcrp.optim(xcrp,fast.only=TRUE)
# print(c(ik, k, il, l, t, b[ik,il,t+1,]))
}
}
}
for (ik in seq_along(ks)) {
for (il in seq_along(ls)) {
w[ik,il,] <- crp(x,b[ik,il,,])
}
}
return(list(b,w))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.