Commit 666202da authored by Martin Maechler's avatar Martin Maechler

add "CLA" class and print() method

parent 74098c55
Functions
^CLA-ex.*Rout
.*\+M$
^SP500-ex.Rout\.
......@@ -11,3 +11,5 @@ export(CLA
, findSig, findMu
, muSigmaGarch
)
S3method(print, CLA)
#### Started from "Version 8" (Ver8.R)
## (for now with minimal changes ( <==> "diff"able)
# Initialize the weight -- Find first free weight
initAlgo <- function(mu, lB, uB ){
#New-ordered return, lB, uB with decreasing return
## Initialize the weight -- Find first free weight
initAlgo <- function(mu, lB, uB ) {
## New-ordered return, lB, uB with decreasing return
w <- c()
index.new <- order(mu,decreasing = TRUE) # new order with decreasing return
lB.new <- lB[index.new]
uB.new <- uB[index.new]
# free weight - starting solution
## free weight - starting solution
i.new <- 0
w.new <- lB.new # initialy
w.new <- lB.new # initially
while(sum(w.new) < 1) {
i.new <- i.new + 1
w.new[i.new] <- uB.new[i.new]
i.new <- i.new + 1
w.new[i.new] <- uB.new[i.new]
}
w.new[i.new] <- 1 - sum(w.new[-i.new])
w[index.new] <- w.new #back to original order
w[index.new] <- w.new # back to original order
i <- index.new[i.new]
## return the index of first free asset and vector w :
list(index = i, weights = w)
}
}
# getMatrices -----------------------------------------------------------------
getMatrices <- function(mu, covar, w, f){
# Slice covarF,covarFB,covarB,muF,muB,wF,wB
## getMatrices -----------------------------------------------------------------
getMatrices <- function(mu, covar, w, f) {
## Slice covarF,covarFB,covarB,muF,muB,wF,wB
covarF <- covar[f,f]
muF <- mu[f]
b <- (seq_along(mu))[-f]
b <- seq_along(mu)[-f]
covarFB <- covar[f,b]
wB <- w[b]
list(covarF = covarF, covarFB = covarFB, muF = muF, wB = wB)
}
}
computeInv <- function(get){
computeInv <- function(get) {
solve(get$covarF, cbind(1, get$muF, get$covarFB %*% get$wB, deparse.level = 0L))
}
}
# computeW -----------------------------------------------------------------
computeW <- function(lam, inv, wB){
# w2 <- inv[,1]; w3 <- inv[,2]; w1 <- inv[,3]
## computeW -----------------------------------------------------------------
computeW <- function(lam, inv, wB) {
## w2 <- inv[,1]; w3 <- inv[,2]; w1 <- inv[,3]
inv.s <- colSums(inv) # g1 <- inv.s[2]; g2 <- inv.s[1]; g4 <- inv.s[3]
#1) compute gamma
## 1) compute gamma
g <- (-lam * inv.s[2] + (1- sum(wB) + inv.s[3]))/inv.s[1]
#2) compute free weights
## 2) compute free weights
list(wF = - inv[,3] + g * inv[,1] + lam * inv[,2], gamma = g)
}
}
# computeLambda --------------------------------------------------------------
computeLambda <- function(wB, inv, i, bi.input){
## computeLambda --------------------------------------------------------------
computeLambda <- function(wB, inv, i, bi.input) {
inv.s <- colSums(inv)
# c1 <- inv.s[1]; l2 <- inv.s[3]; c2i <- inv[i,2];
# c3 <- inv.s[2]; c4i <- inv[i, 1]; l1 <- sum(wB)
## c1 <- inv.s[1]; l2 <- inv.s[3]; c2i <- inv[i,2];
## c3 <- inv.s[2]; c4i <- inv[i, 1]; l1 <- sum(wB)
c1 <- inv.s[1]
if(length(bi.input)==1){ # 1.bound to free
c4i <- inv[i, 1]
Ci <- - c1 * inv[i, 2] + inv.s[2] * c4i
if(Ci == 0) 0
((1- sum(wB) + inv.s[3])* c4i- c1 * (bi.input + inv[i, 3]))/Ci # return lambda
} else { # 2.free to bound
c4i <- inv[, 1]
Ci <- - c1 * inv[, 2] + inv.s[2] * c4i
bi <- bi.input[i, 1] # bi.lB
bi[Ci > 0] <- bi.input[i[Ci > 0], 2] # bi.uB
bi[Ci == 0] <- 0
list(lambda = ((1- sum(wB) + inv.s[3]) * c4i- c1 *(bi + inv[, 3]))/Ci,
bi = bi)
# return lambda and boundary
if(length(bi.input) == 1L) { # 1.bound to free
c4i <- inv[i, 1]
Ci <- - c1 * inv[i, 2] + inv.s[2] * c4i
## return lambda :
if(Ci == 0)
0
else
((1- sum(wB) + inv.s[3])* c4i- c1 * (bi.input + inv[i, 3]))/Ci
}
else { # 2.free to bound
c4i <- inv[, 1]
Ci <- - c1 * inv[, 2] + inv.s[2] * c4i
bi <- bi.input[i, 1] # bi.lB
bi[Ci > 0] <- bi.input[i[Ci > 0], 2] # bi.uB
bi[Ci == 0] <- 0
## return lambda and boundary :
list(lambda = ((1- sum(wB) + inv.s[3]) * c4i- c1 *(bi + inv[, 3]))/Ci,
bi = bi)
}
}
}
MS <- function(weights_set, mu, covar){
MS <- function(weights_set, mu, covar) {
Sig2 <- colSums(weights_set *(covar %*% weights_set) )
cbind(Sig = sqrt(Sig2), Mu = as.vector(t(weights_set) %*% mu))
}
}
CLA <- function(mu, covar, lB, uB, tol.lambda = 1e-7, give.MS = TRUE, keep.names = TRUE) {
CLA <- function(mu, covar, lB, uB, tol.lambda = 1e-7,
give.MS = TRUE, keep.names = TRUE) {
## minimal argument checks
cl <- match.call()
n <- length(mu)
if(length(lB) == 1) lB <- rep.int(lB, n)
if(length(uB) == 1) uB <- rep.int(uB, n)
stopifnot(is.numeric(mu), is.matrix(covar), identical(dim(covar), c(n,n)),
is.numeric(lB), length(lB) == n,
is.numeric(uB), length(uB) == n, lB <= uB)# and in [0,1]
# Compute the turning points, free sets and weights
## Compute the turning points, free sets and weights
ans <- initAlgo(mu, lB, uB)
f <- ans$index
w <- ans$weights
......@@ -95,9 +99,9 @@
free_indices <- list(f)
lam <- 1 # set non-zero lam
while ( lam > 0 && length(f) < length(mu)) {
# 1) case a): Bound one free weight F -> B
## 1) case a): Bound one free weight F -> B
l_in <- 0
if(length(f) > 1 ){
if(length(f) > 1 ) {
compl <- computeLambda(wB = w[-f], inv = inv, # inv from last step k (k > 1)
i = f, bi.input = cbind(lB, uB))
lam_in <- compl$lambda
......@@ -108,9 +112,9 @@
l_in <- lam_in[k]
}
# 2) case b): Free one bounded weight B -> F
## 2) case b): Free one bounded weight B -> F
b <- seq_along(mu)[-f]
inv_list <- lapply(b, function(bi){
inv_list <- lapply(b, function(bi) {
get_i <- getMatrices(mu, covar, w, c(f,bi))
computeInv(get_i)
})
......@@ -131,10 +135,10 @@
l_out <- lam_out[k]
inv_out <- inv_list[[k]]
# 3) decide lambda
## 3) decide lambda
lam <- max(l_in, l_out, 0)
if(lam > 0) { # remove i_in from f; or add i_out into f
if(l_in > l_out ){
if(l_in > l_out ) {
w[i_in] <- bi_in # set value at the correct boundary
f <- f[f != i_in]
getM <- getMatrices(mu, covar, w, f)
......@@ -146,9 +150,9 @@
}
compW <- computeW(lam, inv = inv, wB = w[-f])
}
else{ #4) if max(l_in, l_out) < 0, "stop" when at the min var solution!
else{ # 4) if max(l_in, l_out) < 0, "stop" when at the min var solution!
compW <- computeW(lam = lam, inv = inv, wB = w[-f])
# muF = 0 not necessary, get1 replaced by getM (ie getM from previous step)
## muF = 0 not necessary, get1 replaced by getM (ie getM from previous step)
}
wF <- compW$wF
......@@ -159,12 +163,37 @@
weights_set <- cbind(weights_set, w, deparse.level = 0L) # store solution
gammas <- c(gammas, g)
free_indices <- c(free_indices, list(sort(f)))
} #end While
}# end While
if(keep.names) rownames(weights_set) <- names(mu)
list(weights_set = weights_set,
free_indices = free_indices,
gammas = gammas, lambdas = lambdas,
MS_weights = if(give.MS) MS(weights_set = weights_set, mu = mu, covar = covar))
}
#### FIXME: This change also needs changes in the ../tests/SP500-ex.R !
structure(class = "CLA",
list(weights_set = weights_set,
free_indices = free_indices,
gammas = gammas, lambdas = lambdas,
MS_weights = if(give.MS)
MS(weights_set = weights_set, mu = mu, covar = covar),
call = cl))
}
## print method
print.CLA <- function(x, ...) {
cat("Call:", paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n", sep = "")
wts <- x$weights_set
n <- nrow(wts)
nT <- ncol(wts)
cat(gettextf("Critical Line Algorithm for n = %d assets, resulting in\n",
n),
gettextf("%d turning points\n", nT))
## For now; we can do better later:
cat("Overview of result parts:\n")
utils::str(x)
## TODO better, e.g., summarizing the number "active assets"
## and those with non-0 weights -- only if lower bounds were (mostly) 0
invisible(x)
}
## TODO: plot method -- efficient frontier
......@@ -54,10 +54,12 @@ set.seed(47)
iS <- sample.int(length(muS.sp500$mu), 24)
CLsp.24 <- CLA(muS.sp500$mu[iS], muS.sp500$covar[iS, iS], lB=0, uB=1/10)
CLsp.24 # using the print() method for class "CLA"
if(require(Matrix)) { ## visualize how weights change "along turning points"
image(Matrix(CLsp.24$weights_set, sparse=TRUE),
xlab = "turning point", ylab = "asset number")
show(image(Matrix(CLsp.24$weights_set, sparse=TRUE),
main = "CLA(muS.sp500 <random_sample(size=24)>) $ weights_set",
xlab = "turning point", ylab = "asset number"))
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment