MPT

La Théorie Moderne du Portefeuille (MPT) de H. Markowitz en quelques lignes de code.


rm(list=ls())

# -------------------------------------------------------------------------- #
# yHist
# Pour télécharger des données sur ichart.yahoo.com
# -------------------------------------------------------------------------- #

yHist = function(sym, from, to = Sys.Date(), freq = "d") {
      f <- as.Date(from[1])
      t <- as.Date(to[1])
      q <- tolower(substr(freq[1], 1, 1))
      if(! q %in% c("d", "w", "m")) stop("! freq %in% c('d', 'w', 'm')")
      cn <- c("Sym", "Date", "Adj.Close")
      if(length(sym) > 1) {
            k <- lapply(sym, yHist, from = from, to = to, freq = freq)
            res <- do.call(rbind, k)
      } else {
            url <- paste(c(
                 "http://ichart.yahoo.com/table.csv?",
                  "s=", sym,
                  "&a=", as.integer(format(f, "%m")) - 1,
                  "&b=", as.integer(format(f, "%d")),
                  "&c=", as.integer(format(f, "%Y")),
                  "&d=", as.integer(format(t, "%m")) - 1,
                  "&e=", as.integer(format(t, "%d")),
                  "&f=", as.integer(format(t, "%Y")),
                  "&g=", q,
                  "&ignore=.csv"), collapse = "")
            ans <- try(read.csv(url), silent = TRUE)
            if(class(ans) == "try-error") {
                  res <- NULL
            } else {
                  ans <- ans[order(ans$Date), ]
                  res <- subset(cbind(Sym = sym, ans), select = cn)
                  rownames(res) <- NULL
            }
      }
      return(res)
}

# -------------------------------------------------------------------------- #
# dvar
# Calcule des % de variation.
# -------------------------------------------------------------------------- #

dvar = function(x) {
 if( is.matrix(x) ) {
  res <- apply(x, 2, dvar)
  dimnames(res) <- list(rownames(x)[-1], colnames(x))
 } else {
  res <- x[-1]/x[-length(x)]-1
  names(res) <- names(x)[-1]
 }
 return(res)
}

# ========================================================================== #
# SCRIPT
# ========================================================================== #

# L'indice DJIA (j'ai juste remplacé GS par AAPL).

Sym <- c("AAPL", "AXP", "BA", "CAT", "CSCO", "CVX", "DD", "DIS", "GE", "HD",
 "IBM", "INTC", "JNJ", "JPM", "KO","MCD", "MMM", "MRK", "MSFT", "NKE",
 "PFE", "PG", "T", "TRV", "UNH", "UTX", "V", "VZ", "WMT", "XOM")

# Téléchargement

Dta <- yHist(c(Sym, "^GSPC"), "2013-11-20", "2014-11-20", "d")

# On transforme tout ça en matrice de prix.

P <- with(Dta, tapply(Adj.Close, list(Date, Sym), mean))

# On sauvegarde les dates pour usage ultérieur.

Dates <- as.Date(rownames(P))

# On isole le S&P 500 du reste de nos données.

Ix <- P[, "^GSPC", drop = FALSE]
P <- P[, 1:(ncol(P)-1)]

# Un taux sans risque bricolé rapidement...

Rf <- as.matrix(1.0012^(as.numeric(diff(Dates))/365)-1)

# Matrices de variations quotidiennes

X <- dvar(P)
Xm <- dvar(Ix)

# Portefeuille équipondéré

w <- matrix(1/30, 30, 1)
rownames(w) <- Sym

# Qu'est-ce que ça donne ?
po <- X %*% w

y <- c(1, cumprod(1 + po))
plot(Dates, y, type = "l")

# Voyons la volatilité de ce portefeuille.
sd(po) * sqrt(260)

# ========================================================================== #
# MINIMUM VARIANCE
# ========================================================================== #

# On charge *quadprog*. Si vous ne l'avez pas : install.packages("quadprod")
library("quadprog")

# Matrice de variance-covariance.
S <- cov(X)

# Tout le monde à 0 !
R <- matrix(0, 1, 30)

# Nous voulons que la somme des poids soit égale à 1 ET que tous les poids
# soient supérieurs à 0.
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

# Attention !
res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# 21% de McDo
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Ce portefeuille est imbattable
y <- c(1, cumprod(1 + po))
plot(Dates, y, type = "l")

# ========================================================================== #
# INDICIEL
# ========================================================================== #

Xs <- X - drop(Xm)

# Matrice de variance-covariance.
S <- cov(Xs)

# Tout le monde à 0 !
R <- matrix(0, 1, 30)

# Nous voulons que la somme des poids soit égale à 1 ET que tous les poids
# soient supérieurs à 0.
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

# Attention !
res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

po <- X %*% w
sd(po - Xm) * sqrt(260)

# Visualisons
y <- c(1, cumprod(1 + po))
ym <- c(1, cumprod(1 + Xm))
plot(Dates, y, type = "l")
lines(Dates, ym, lty = "dashed")

# ========================================================================== #
# LA FRONTIERE EFFICIENTE
# ========================================================================== #

R <- colMeans(X)
S <- cov(X)
A <- cbind(1, diag(1, 30, 30))
B <- matrix(c(1, rep(0, 30)))

res <- solve.QP(S, R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# LOL ! 100% Apple !
round(w*100, 2)

# Comment faire baisser cette volatilité ?
# Facile : on va utiliser le paramètre q

q <- .5
res <- solve.QP(S, q*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# Intel rentre dans notre portefeuille.
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Réduisons encore q :

q <- .1
res <- solve.QP(S, q*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

# UnitedHealth rentre dans notre portefeuille.
round(w*100, 2)

po <- X %*% w
sd(po) * sqrt(260)

# Voyons ce que ça donne en faisant varier q entre 0 et .1 :
q <- seq(0, .1, len = 100)

W1 <- matrix(NA, 100, 30)
colnames(W1) <- Sym

Z1 <- matrix(NA, 100, 2)
colnames(Z1) <- c("mean", "sd")

for(i in 1:length(q)) {
 ans <- solve.QP(S, q[i]*R, A, B, 1)
 W1[i, ] <- ans$solution
 po <- X %*% matrix(ans$solution, 30, 1)
 Z1[i, ] <- c(mean(po), sd(po))
}

# La frontière efficiente !
plot(Z1[, "sd"], Z1[, "mean"], type = "b")

# ========================================================================== #
# ACTIF SANS RISQUE
# ========================================================================== #

XX <- cbind(X, Rf)
R <- colMeans(XX)
S <- cov(XX)
A <- cbind(1, diag(1, 31, 31))                 # Attention : 31 !
B <- matrix(c(1, rep(0, 31)))                  # Attention : 31 !

q <- seq(0, .1, len = 100)

W2 <- matrix(NA, 100, 31)                      # Attention : 31 !
colnames(W2) <- Sym

Z2 <- matrix(NA, 100, 2)
colnames(Z2) <- c("mean", "sd")

for(i in 1:length(q)) {
 ans <- solve.QP(S, q[i]*R, A, B, 1)
 W2[i, ] <- ans$solution
 po <- XX %*% matrix(ans$solution, 31, 1) # Attention : 31 !
 Z2[i, ] <- c(mean(po), sd(po))
}

# Nouvelle frontière efficiente : la Capital Allocation Line !
plot(Z2[, "sd"], Z2[, "mean"], type = "b")

# ========================================================================== #
# LONG SHORT !
# ========================================================================== #

R <- colMeans(X)
S <- cov(X)
A <- cbind(1, diag(1, 30, 30), diag(-1, 30, 30))
B <- matrix(c(0, rep(-.2, 30), rep(-.2, 30)))

res <- solve.QP(S, .01*R, A, B, 1)

w <- matrix(res$solution, 30, 1)
rownames(w) <- Sym

po <- X %*% w
sd(po) * sqrt(260)

Aucun commentaire:

Enregistrer un commentaire