8 Portfolioptimierung
Wir betrachten im Folgenden, wie mittels R die quadratische (Mean-Variance) und die lineare (Mean-Expected-Shortfall) Optimierung durchgeführt werden kann. In beiden Fällen wird ein entsprechender Algorithmus zur Lösung derartiger Probleme verwendet. Wichitg ist es bei der Implementierung darauf zu achten, wie das Optimierungsproblem übergeben werden muss, so dass es vom Algorithmus gelöst werden kann. Für beide Fälle beziehen wir als Beispiel Daten das Dax. In der Realität ist es jedoch nicht sinnvoll lediglich in einem Markt zu investieren, da auf internationaler und vor allem interkontinentaler Ebene hohes Diversifikationspotential besteht.
rm(list = ls())
library(tidyquant)
library(data.table)
library(plotly)
library(quadprog)
library(FRAPO)
library(Rglpk)
library(VineCopula)
library(slam)
###########################################################################################
#Get the DAX data #
###########################################################################################
#yahoo finance tickers
yh_tickers <- read.csv('/Users/ralfkellner/Library/Mobile Documents/com~apple~CloudDocs/Lehre Saarbrücken/QuantMeth/Codes/Markdown/Data/dax_yahoo_tickers.csv')
tickers <- as.character(yh_tickers$Ticker)
#data
our_data <- tq_get(tickers, from = "2010-02-01", to = "2020-02-03")
#transform to data table
df <- data.table(our_data)
#transform to wide format with closing values
dt_close <- dcast(df, date ~ symbol , value.var = "close")
head(dt_close)## date 1COV.DE ADS.DE ALV.DE BAS.DE BAYN.DE BEI.DE BMW.DE CON.DE DAI.DE
## 1: 2010-02-01 NA 37.065 80.49 41.730 48.8210 42.275 31.065 41.320 33.570
## 2: 2010-02-02 NA 36.975 81.19 42.000 49.2048 41.965 31.175 41.280 33.765
## 3: 2010-02-03 NA 36.280 80.41 42.130 48.4126 42.220 31.225 41.100 34.340
## 4: 2010-02-04 NA 35.560 79.29 40.255 47.5810 41.795 30.335 39.030 33.100
## 5: 2010-02-05 NA 35.175 77.61 39.430 46.0753 41.275 29.920 37.765 32.320
## 6: 2010-02-08 NA 35.005 78.11 40.020 46.7986 41.910 29.610 36.400 33.310
## DB1.DE DBK.DE DPW.DE DTE.DE EOAN.DE FME.DE FRE.DE HEI.DE HEN3.DE IFX.DE
## 1: 47.965 35.4626 12.810 9.410 24.3161 36.685 14.7067 43.375 36.26 4.208
## 2: 48.950 36.2355 13.145 9.440 24.5744 37.130 14.6700 44.530 36.12 4.343
## 3: 47.240 35.5869 12.925 9.342 24.2844 36.990 14.6767 44.170 35.83 4.220
## 4: 47.095 34.0839 12.415 9.240 23.7905 36.225 14.5717 42.600 35.21 4.110
## 5: 45.860 33.5440 12.140 9.211 23.2014 36.115 14.3833 41.420 35.29 3.959
## 6: 46.105 33.4780 12.210 9.372 23.5186 36.380 14.7050 41.935 35.81 3.958
## LHA.DE LIN.DE MRK.DE MUV2.DE RWE.DE SAP.DE SIE.DE TKA.DE VNA.DE VOW3.DE
## 1: 11.770 80.18 32.065 109.35 64.7643 33.585 63.0354 23.365 NA 60.2575
## 2: 11.980 80.65 32.160 108.75 65.5318 33.795 64.0529 23.785 NA 60.6750
## 3: 11.745 80.89 32.060 107.95 64.5750 34.050 63.4424 23.550 NA 61.1323
## 4: 11.330 79.56 31.755 106.90 63.2195 33.795 61.0973 22.350 NA 60.2377
## 5: 11.045 79.26 31.275 106.25 62.4421 33.385 59.7601 22.035 NA 58.8858
## 6: 11.080 79.74 31.570 107.55 63.3690 32.560 60.6225 22.280 NA 58.2993
## WDI.DE
## 1: 9.208
## 2: 9.570
## 3: 9.173
## 4: 9.048
## 5: 8.820
## 6: 8.756#eliminate those with NAs over time
cols <- names(dt_close)[(apply(dt_close, 2, function(x){sum(is.na(x))}) == 0)]
dt_close <- dt_close[ , cols, with = FALSE]
head(dt_close)## date ADS.DE ALV.DE BAS.DE BAYN.DE BEI.DE BMW.DE CON.DE DAI.DE DB1.DE
## 1: 2010-02-01 37.065 80.49 41.730 48.8210 42.275 31.065 41.320 33.570 47.965
## 2: 2010-02-02 36.975 81.19 42.000 49.2048 41.965 31.175 41.280 33.765 48.950
## 3: 2010-02-03 36.280 80.41 42.130 48.4126 42.220 31.225 41.100 34.340 47.240
## 4: 2010-02-04 35.560 79.29 40.255 47.5810 41.795 30.335 39.030 33.100 47.095
## 5: 2010-02-05 35.175 77.61 39.430 46.0753 41.275 29.920 37.765 32.320 45.860
## 6: 2010-02-08 35.005 78.11 40.020 46.7986 41.910 29.610 36.400 33.310 46.105
## DBK.DE DPW.DE DTE.DE EOAN.DE FME.DE FRE.DE HEI.DE HEN3.DE IFX.DE LHA.DE
## 1: 35.4626 12.810 9.410 24.3161 36.685 14.7067 43.375 36.26 4.208 11.770
## 2: 36.2355 13.145 9.440 24.5744 37.130 14.6700 44.530 36.12 4.343 11.980
## 3: 35.5869 12.925 9.342 24.2844 36.990 14.6767 44.170 35.83 4.220 11.745
## 4: 34.0839 12.415 9.240 23.7905 36.225 14.5717 42.600 35.21 4.110 11.330
## 5: 33.5440 12.140 9.211 23.2014 36.115 14.3833 41.420 35.29 3.959 11.045
## 6: 33.4780 12.210 9.372 23.5186 36.380 14.7050 41.935 35.81 3.958 11.080
## LIN.DE MRK.DE MUV2.DE RWE.DE SAP.DE SIE.DE TKA.DE VOW3.DE WDI.DE
## 1: 80.18 32.065 109.35 64.7643 33.585 63.0354 23.365 60.2575 9.208
## 2: 80.65 32.160 108.75 65.5318 33.795 64.0529 23.785 60.6750 9.570
## 3: 80.89 32.060 107.95 64.5750 34.050 63.4424 23.550 61.1323 9.173
## 4: 79.56 31.755 106.90 63.2195 33.795 61.0973 22.350 60.2377 9.048
## 5: 79.26 31.275 106.25 62.4421 33.385 59.7601 22.035 58.8858 8.820
## 6: 79.74 31.570 107.55 63.3690 32.560 60.6225 22.280 58.2993 8.756#derive returns in %
returns <- apply(log(dt_close[, -c(1)]), 2, diff) * 100
#show correlations in heatmap
corr <- cor(returns)
fig <- plot_ly(x = rownames(corr), y = rownames(corr), z = corr, type = "heatmap")
fig8.1 Quadratische Optimierung
Als Inputdaten bei der Optimierung schätzen wir die Erwartungswerte der Renditen und die Kovarianzmatrix über ihre empirichen Schätzer.
#estimate means
means <- apply(stocks, 2, mean)
#estimate sds
sds <- apply(stocks, 2, sd)
#estimated covariance matrix
cov_est <- cov(stocks)
n_assets <- dim(stocks)[2]Die Optimierung wird mittels der solve.QP Funktion durchgeführt. In dieser wird neben der Kovarianzmatrix als Zielfunktion eine Matrix und ein Vektor mit Nebenbedingungen angegeben, wobei der meq-Parameter der Funktion angibt, wie viele der ersten Nebenbedingungen Gleicheitsbedingungen sind. Bei den übrigen Nebenbedingungen geht die Funktion von Ungleichungen mit \(\geq\) Relation aus. Möchten wir beispielsweise lediglich volle Investition und einen Mindesrendite als Nebenbedingung, so lautet die Matrix und der Vektor.
#constraint matrix
A_mat <- rbind(rep(1, n_assets), #for full investment
means) #for desired portfolio return
#conditions
b_vec <- c(1, 0.025)
#solving quadratic problem
quad_opt <- try(solve.QP(Dmat = 2 * cov_est, dvec = rep(0, n_assets),
Amat = t(A_mat), bvec = b_vec, meq = 2))
quad_opt## $solution
## [1] -0.005890768 -0.159815940 -0.040064640 -0.055556610 0.255169959
## [6] 0.017720191 -0.043403970 -0.012618238 0.079585543 0.025115598
## [11] 0.055758064 0.103222661 0.056010272 0.139912387 0.014012396
## [16] -0.030463805 0.075166343 -0.076075730 0.037107148 0.170960541
## [21] 0.050711494 0.185741384 -0.006829584 0.106026785 0.082627930
## [26] 0.001540095 -0.001842687 -0.023826818
##
## $value
## [1] 0.7431854
##
## $unconstrained.solution
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##
## $iterations
## [1] 3 0
##
## $Lagrangian
## [1] 1.576049 3.587134
##
## $iact
## [1] 1 2Um effiziente Portfolioallokationen für verschiedene Niveaus der Rendite zu bestimmen, führen wir die Optimierung in einer Schleife durch. Zuletzt bestimmen wir die Sharpe-Ratio, wobei wir vereinfachend aufgrund des aktuellen Zinsniveaus \(r_f = 0.00\) verwenden.
#desired expected portfolio return
mu_pf <- seq(0, 0.10, length.out = 50)
#constraint matrix
A_mat <- rbind(rep(1, n_assets), #for full investment
means, #for desired portfolio return
diag(1, n_assets), #set minimum weight to invest
diag(-1, n_assets)) #set maximum weight to invest
#record weights and standard deviation
weights <- matrix(ncol = n_assets)
sd_pf <- c()
min_weight <- -1
max_weight <- 1
for(i in 1:length(mu_pf)){
#constraint vecotr
b_vec <- c(1, mu_pf[i], rep(min_weight, n_assets), rep(-max_weight, n_assets))
#solving quadratic problem
quad_opt <- try(solve.QP(Dmat = 2 * cov_est, dvec = rep(0, n_assets),
Amat = t(A_mat), bvec = b_vec, meq = 2))
#some solutions may not be attainable under constraints
#if this is the case save NA
if(class(quad_opt) == 'try-error'){
weights <- rbind(weights, rep(NA, n_assets))
sd_pf[i] <- NA
}else{
#save results
weights <- rbind(weights, quad_opt$solution)
sd_pf[i] <- quad_opt$value
}
}
p <- plot_ly()%>%
add_markers(x = sd_pf, y = mu_pf, name = 'efficient frontier', marker = list(color = 'grey'))%>%
layout(xaxis = list(title = 'standard deviation', range = c(0, max(sd_pf))),
yaxis = list(title = 'mean'))
r_f <- 0
sharpe <- (mu_pf - r_f) / sd_pf
ind <- match(max(sharpe, na.rm = TRUE), sharpe)
weights[ind,]## [1] 0.143295028 0.134141422 -0.101835182 -0.101658612 0.171340274
## [6] 0.096540545 -0.024689837 -0.170414097 0.115931222 -0.180007548
## [11] 0.111507991 0.050530340 -0.125584095 0.015797299 0.084986219
## [16] -0.059302066 0.053796890 0.014666841 -0.001367922 0.141248452
## [21] 0.155447312 0.260361471 -0.027854546 0.146839003 0.068272252
## [26] -0.095174359 0.059255632 0.0639300728.2 Lineare Optimierung
Ganz analog kann die lineare Optimierung durch die Rglpk_solve_LP Funktion durchgeführt werden. Man sollte hierzu im Kopf haben, dass zwar nur über \(\boldsymbol{\omega}, \psi\) optimiert wird, jedoch auch für \(z_j, ~j = 1, ..., J\) Werte bei der Optimierung mit bestimmt werden. Daher gilt es diese beim Festlegen der Ziel- und Nebenfunktionen mit zu berücksichtigen. Eine detaillierte Erläuterung zur Verwendung der Funktion erfolgt in einem Lehrvideo.
n_assets <- ncol(stocks)
n_sample <- nrow(stocks)
conf_level <- 0.95
#objective function (asset weights, ..., z_1, z_2, ..., z_n, psi)
obj <- c(rep(0, n_assets), rep(1 / ((1 - conf_level) * n_sample), n_sample), 1)
#full investment #simple_triplet_format to save memory space
M1 <- cbind(matrix(rep(1, n_assets), nrow = 1, ncol = n_assets) , simple_triplet_zero_matrix(nrow = 1, ncol = n_sample + 1))
# -f(w_j, x_j) + alpha + z_j >= 0
M2 <- cbind(stocks, simple_triplet_diag_matrix(1, n_sample), matrix(1, nrow = n_sample) )
#z_j >= 0
M3 <- cbind(simple_triplet_zero_matrix(ncol = n_assets, nrow = n_sample), simple_triplet_diag_matrix(1, n_sample), matrix(0, nrow = n_sample))
#desired portfolio return
M4 <- cbind(matrix(apply(stocks, 2, mean), nrow = 1, ncol = n_assets), simple_triplet_zero_matrix(nrow = 1, ncol = n_sample + 1))
#min weights
M5 <- cbind(simple_triplet_diag_matrix(1, nrow = n_assets), simple_triplet_zero_matrix(nrow = n_assets, ncol = n_sample + 1))
#max weights
M6 <- cbind(simple_triplet_diag_matrix(-1, nrow = n_assets), simple_triplet_zero_matrix(nrow = n_assets, ncol = n_sample + 1))
#constraints matrix
Amat <- rbind(M1, M2, M3, M4, M5, M6)
#directions of constraints
dirVec <- c("==", rep(">=", 2*n_sample), "==", rep(">=", n_assets), rep(">=", n_assets))
#bounds for optimization parameters
bounds <- list(lower = list(ind = c(1:n_assets,(n_assets+1):(n_assets+n_sample)), val = c(rep(0,n_assets),rep(-10,n_sample))),
upper = list(ind = c(1:n_assets,(n_assets+1):(n_assets+n_sample)), val = c(rep(1,n_assets),rep(10,n_sample))))
#record weights, CVaR and corresponding standard deviation
weights_cvar <- matrix(ncol = n_assets)
CVaR <- c()
sd_CVaRpf <- c()
#desired expected portfolio return
mu_pf <- seq(0, 0.10, length.out = 50)
#minimum weights
min_weight <- -1
#maximum weights #max weight constraint is introduced as -w_j >= -max_weight <==> max_weight >= w_j
max_weight <- 1
for(i in 1:length(mu_pf)){
#constraints vector
b_vec <- c(1, rep(0, 2 * n_sample), mu_pf[i], rep(min_weight, n_assets), rep(-max_weight, n_assets))
#linear optimization routine
results <- try(Rglpk_solve_LP(obj = obj, mat = Amat, dir = dirVec, rhs = b_vec, bounds = bounds, max = FALSE))
#some solutions may not be attainable under constraints
#if this is the case save NA
if(class(results) == 'try-error'){
CVaR[i] <- NA
sd_CVaRpf[i] <- NA
weights_cvar <- rbind(weights_cvar, rep(NA, n_assets))
}else{
#results
r_opt <- as.vector(stocks %*% matrix(results$solution[1:n_assets], nrow = n_assets))
#VaR <- tail(results$solution, 1)
CVaR[i] <- results$optimum
sd_CVaRpf[i] <- sd(r_opt)
weights_cvar <- rbind(weights_cvar, results$solution[1:n_assets])
}
#print(i)
}
p <- plot_ly()%>%
add_markers(x = CVaR, y = mu_pf, name = 'efficient frontier', marker = list(color = 'grey'))%>%
layout(xaxis = list(title = 'expected shortfall', range = c(0, max(CVaR))),
yaxis = list(title = 'mean'))
r_f <- 0
sharpe_CVaR <- (mu_pf - r_f) / CVaR
ind <- match(max(sharpe_CVaR, na.rm = TRUE), sharpe_CVaR)
weights_cvar[ind,]## [1] 0.48124322 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [7] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [13] 0.00000000 0.00000000 0.10100899 0.00000000 0.00000000 0.00000000
## [19] 0.00000000 0.05579814 0.07679953 0.00000000 0.00000000 0.09984039
## [25] 0.00000000 0.00000000 0.00000000 0.18530972