In AxelCode-R/Master-Thesis: Does not matter.
Deterministic Optimization Methods for Quadratic Programming Problems {#analyticalsolver}
\chaptermark{Deterministic Optimization Methods}
In this chapter, the advantages and disadvantages of deterministic optimization methods for quadratic programming problems are discussed. For simplicity, this thesis refers to the implementation of a deterministic optimization method as a solver. It is beyond the scope of this thesis to explain the underlying mathematical principles of how a solver handles quadratic problems; only the applications and analysis are discussed. The main reason for dealing with solvers for quadratic programming problems is to use them as a benchmark for the PSO.
Quadratic Programming (QP)
A quadratic program is a minimization problem of a function that returns a scalar and consists of a quadratic term and a linear term that depend on the variable of interest. In addition, the problem may be constrained by several linear inequalities that bound the solution. The general formulation used is to find $x$ that minimizes the following problem:
$$
\min\limits_{x} \ \frac{1}{2} \cdot x^T \times D \times x - d^T \times x
$$
and is valid under the linear constraints:
$$
A^T \times x \geq b_0.
$$
Some other sources notate the problem with different signs or coefficients, all of which are interchangeable with the above problem. In addition, the above problem has the same notation as in the R package quadprog, which reduces the substitution overhead. All modern programming languages have many solvers for quadratic problems. They differ mainly in the computation time for certain problems and the requirements. Some commercial QP solvers additionally accept more complex constraints, such as absolute (e.g., $|A^T \times x| \geq a_0$) or mixed-integer (e.g., $x \in \mathbb{N}$). Especially the mixed-integer constraint problems lead to a huge increase in memory requirements.
QP Solver from quadprog
The common free QP solver used in R comes from the package quadprog, which consists of a single function called solve.QP. Its implementation routine is the dual method of Goldfarb and Idnani published in [@GoId1982] and [@GoId1983]. It uses the above QP with the condition that $D$ has to be a symmetric positive definite matrix. This means that $D\in \mathbb{R}^{N \times N}$ and $x^T D x > 0 \ \forall \ x \in \mathbb{R}^N \setminus {\vec{0}}$, which is equivalent to all eigenvalues being greater than zero. In most cases this is not achieved by estimating the covariance matrix $\sum$, but it is possible to find the nearest positive definite matrix of $\textstyle\sum$ using the function nearPD() from the matrix R package. The error encountered often does not exceed a percent change in elements greater than $10^{-15} \%$, which is negligible for the context of this work. The function solve.QP for an $N$ dimensional vector of interest, has the following arguments, which are also found in the above formulation of a QP:
Dmat: Symmetric positive definite matrix $D \in \mathbb{R}^{N \times N}$ of the quadratic termdvec: Vector $d \in \mathbb{R}^{N}$ of the linear termAmat: Constraint matrix $A$bvec: Constraint vector $b_0$meq = 1: means that the first meq columns of $A$ are treated as an equality constraint
The return of solve.QP is a list and contains, among others, the following attributes of interest:
solution: Vector containing the solution $x$ of the quadratic programming problem (e.g. portfolio weights)value: Scalar, the value of the quadratic function at the solution
Solving MVP Problem with solve.QP {#exampleanalyticalmvp}
This section provides insights into the effects of diversification and the use of solve.QP by creating ten different efficient frontiers from a pool of ten assets. Each efficient frontier $i \in {1, 2, \cdots, 10}$ consists of $N_i = i$ assets and is created by adding the asset with the next smallest variance first. After loading the returns for ten of the largest stocks in the U.S. market, the variance is calculated to rank all columns in ascending order of variance, as shown in the code below:
returns_raw <- buffer(
get_yf(
tickers = c("IBM", "GOOG", "AAPL", "MSFT", "AMZN",
"NVDA", "JPM", "META", "V", "WMT"),
from = "2018-01-01",
to = "2019-12-31"
)$returns,
"AS_10_assets"
)
# re-arrange: low var first
vars <- sapply(returns_raw, var)
returns_raw <- returns_raw[, order(vars, decreasing = F)]
The next step is to create a function mvp() that has the arguments return and lambda. It computes the expected returns mu and the estimated positive definite covariance cov. It then solves an MVP problem with constraints $\textstyle\sum w_i = 1$ and $w_i \geq 0$, which yields the key features mu, var and composition of the portfolio.
mvp <- function(returns, lambda){
tc <- tryCatch({
mu <- ret_to_geomeanret(returns)
cov <- as.matrix(nearPD(cov_(returns, mu))$mat)
mat <- list(
Dmat = lambda * cov,
dvec = (1-lambda) * mu,
Amat = t(rbind(
rep(1, ncol(returns)), # sum up to 1
diag(
1, nrow=ncol(returns), ncol=ncol(returns)
) # long only
)),
bvec = c(
1, # sum up to 1
rep(0, ncol(returns)) # long only
),
meq = 1
)
qp <- solve.QP(
Dmat = mat$Dmat, dvec = mat$dvec,
Amat = mat$Amat, bvec = mat$bvec, meq = mat$meq
)
res <- list(
"mu" = mu %*% qp$solution,
"var" = t(qp$solution) %*% cov %*% qp$solution,
"composition" = setNames(qp$solution, colnames(returns))
)
TRUE
}, error = function(e){FALSE})
if(tc){
return(res)
}else{
return(list(
"mu" = NA,
"var" = NA,
"composition" = NA
))
}
}
Each $\lambda \in {0.01, 0.02, \cdots, 1}$ and each combination of ascending number of assets results in a portfolio that can be created with two for loops.
df <- data.frame(
"index"=1,
"var"=as.numeric(var(returns_raw[, 1])),
"return" = as.numeric(ret_to_geomeanret(returns_raw[, 1])),
row.names=NULL
)
for(i in 2:ncol(returns_raw)){
returns <- returns_raw[, 1:i]
for(lambda in seq(0.01, 1, 0.01)){
res <- mvp(returns, lambda)
df <- rbind(
df,
data.frame("index"=i, "var"=res$var, "return" = res$mu)
)
}
}
The result is filtered and names are added to represent the number of assets. Now the $\mu$-$\sigma$ diagram can be created:
df <- df %>%
filter(!is.na(return)) %>%
distinct() %>%
mutate(name = paste0("n_", index)) %>%
arrange(name) %>%
mutate(name = factor(name, levels=paste0("n_", ncol(returns_raw):1)))
max_show_sd <- df %>%
group_by(index) %>%
summarise(max_x = max(var)) %>%
pull(max_x) %>%
mean() %>%
sqrt()
plot_ly(
data = df[df$index!=1,],
x=~sqrt(var),
y=~return,
name=~name,
mode="lines",
type = 'scatter',
color = ~name,
colors = c("green", "red")
) %>%
add_trace(
data=df[df$index==1,],
x=~sqrt(var),
y=~return,
showlegend=T,
marker=list(color="red"),
mode="markers",
name="n_1") %>%
layout(
xaxis=list(range=c(sqrt(min(df$var))*0.9, max_show_sd), title="sigma"),
yaxis=list(range=c(min(df$return)*0.9, (max(df$return)+mean(df$return))*0.5), title="mu")) %>%
html_save()
It can be seen that each added asset leads to a minimum variance portfolio with a smaller or equal standard deviation. Despite starting with the asset that has the smallest standard deviation of r sqrt(min(vars)). This is the effect of diversification mentioned by Markowitz.
Solving ITP-MSTE with solve.QP {#exampleitpsolveqp}
This example analyzes how many assets are needed to minimize the mean square error between the replication and historical returns of the SP500TR from 2018-01-01 to 2019-12-31. The constraints are set to be long only and the weights should sum to one. To gradually reduce the number of assets, the five assets with the lowest weights are discarded, and the remaining assets serve as the new asset pool for the next replication until only five assets remain. First, the required data can be downloaded from the R/ directory using existing functions. The function get_spx_composition() uses web scraping to read the components of wikipedia and converts them into monthly compositions of the SP500TR. The pool is formed from all assets present in the last month of the time frame, reduced by assets with missing values. The code below loads the returns of all assets in the pool and the SP500TR:
from <- "2018-01-01"
to <- "2019-12-31"
spx_composition <- buffer(
get_spx_composition(),
"AS_spx_composition"
)
pool_returns_raw <- buffer(
get_yf(
tickers = spx_composition %>%
filter(Date<=to) %>%
filter(Date==max(Date)) %>%
pull(Ticker),
from = from,
to = to
)$returns,
"AS_sp500_assets"
)
pool_returns_raw <-
pool_returns_raw[, colSums(is.na(pool_returns_raw))==0]
bm_returns <- buffer(
get_yf(tickers = "^SP500TR", from = from, to = to)$returns,
"AS_sp500tr"
) %>% setNames(., "SP500TR")
The required data is now available and the function for the ITP-MSTE can be created. It requires pool_returns with variable number of columns and the single-column matrix bm_returns.
itp <- function(pool_returns, bm_returns){
mat <- list(
Dmat = t(pool_returns) %*% pool_returns,
dvec = t(pool_returns) %*% bm_returns,
Amat = t(rbind(
rep(1, ncol(pool_returns)), # sum up to 1
diag(1,
nrow=ncol(pool_returns),
ncol=ncol(pool_returns)) # long only
)),
bvec = c(
1, # sum up to 1
rep(0, ncol(pool_returns)) # long only
),
meq = 1
)
qp <- solve.QP(
Dmat = mat$Dmat, dvec = mat$dvec,
Amat = mat$Amat, bvec = mat$bvec, meq = mat$meq
)
res <- list(
"var" = as.numeric(
var(pool_returns %*% qp$solution - bm_returns)),
"solution" = setNames(qp$solution, colnames(pool_returns))
)
}
The successive removal of assets can begin and the results are stored in res.
res <- NULL
asset_pool <- NULL
n_assets <- rev(seq(5, ncol(pool_returns_raw), 5))
for(i in n_assets){
temp <- if(i==max(n_assets)){
itp(pool_returns_raw, bm_returns)
}else{
asset_pool <- names(sort(temp$solution, decreasing = T)[1:i])
itp(
pool_returns_raw[, asset_pool],
bm_returns
)
}
res <- rbind(
res,
data.frame("N"=i, "var"=temp$var, "sd"=sqrt(temp$var), row.names = NULL)
)
# for later analysis
if(length(asset_pool)==100){
assets_pool_100 <- asset_pool
save(assets_pool_100, file="data/assets_pool_100.rdata")
}
if(length(asset_pool)==50){
assets_pool_50 <- asset_pool
save(assets_pool_50, file="data/assets_pool_50.rdata")
}
}
The following chart illustrates the increase in standard deviation to track the historical performance of the SP500TR as the number of assets $N$ in the asset pool decreases:
plot_ly(data=res, x=~N, y=~sd, mode="lines", type = 'scatter') %>%
layout(yaxis=list(range=c(0, mean(max(res$sd),mean(res$sd)) ))) %>%
html_save()
It can be seen that the standard deviation stagnates at about $N=200$. This leads to the conclusion that a sparse replication with two hundred assets is sufficient in this particular case to track the historical performance of the SP500TR over this period.
AxelCode-R/Master-Thesis documentation built on Feb. 25, 2023, 7:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Deterministic Optimization Methods for Quadratic Programming Problems {#analyticalsolver}
\chaptermark{Deterministic Optimization Methods} In this chapter, the advantages and disadvantages of deterministic optimization methods for quadratic programming problems are discussed. For simplicity, this thesis refers to the implementation of a deterministic optimization method as a solver. It is beyond the scope of this thesis to explain the underlying mathematical principles of how a solver handles quadratic problems; only the applications and analysis are discussed. The main reason for dealing with solvers for quadratic programming problems is to use them as a benchmark for the PSO.
Quadratic Programming (QP)
A quadratic program is a minimization problem of a function that returns a scalar and consists of a quadratic term and a linear term that depend on the variable of interest. In addition, the problem may be constrained by several linear inequalities that bound the solution. The general formulation used is to find $x$ that minimizes the following problem:
$$ \min\limits_{x} \ \frac{1}{2} \cdot x^T \times D \times x - d^T \times x $$
and is valid under the linear constraints:
$$ A^T \times x \geq b_0. $$
Some other sources notate the problem with different signs or coefficients, all of which are interchangeable with the above problem. In addition, the above problem has the same notation as in the R package quadprog, which reduces the substitution overhead. All modern programming languages have many solvers for quadratic problems. They differ mainly in the computation time for certain problems and the requirements. Some commercial QP solvers additionally accept more complex constraints, such as absolute (e.g., $|A^T \times x| \geq a_0$) or mixed-integer (e.g., $x \in \mathbb{N}$). Especially the mixed-integer constraint problems lead to a huge increase in memory requirements.
QP Solver from quadprog
The common free QP solver used in R comes from the package quadprog, which consists of a single function called solve.QP. Its implementation routine is the dual method of Goldfarb and Idnani published in [@GoId1982] and [@GoId1983]. It uses the above QP with the condition that $D$ has to be a symmetric positive definite matrix. This means that $D\in \mathbb{R}^{N \times N}$ and $x^T D x > 0 \ \forall \ x \in \mathbb{R}^N \setminus {\vec{0}}$, which is equivalent to all eigenvalues being greater than zero. In most cases this is not achieved by estimating the covariance matrix $\sum$, but it is possible to find the nearest positive definite matrix of $\textstyle\sum$ using the function nearPD() from the matrix R package. The error encountered often does not exceed a percent change in elements greater than $10^{-15} \%$, which is negligible for the context of this work. The function solve.QP for an $N$ dimensional vector of interest, has the following arguments, which are also found in the above formulation of a QP:
Dmat: Symmetric positive definite matrix $D \in \mathbb{R}^{N \times N}$ of the quadratic termdvec: Vector $d \in \mathbb{R}^{N}$ of the linear termAmat: Constraint matrix $A$bvec: Constraint vector $b_0$meq = 1: means that the firstmeqcolumns of $A$ are treated as an equality constraint
The return of solve.QP is a list and contains, among others, the following attributes of interest:
solution: Vector containing the solution $x$ of the quadratic programming problem (e.g. portfolio weights)value: Scalar, the value of the quadratic function at the solution
Solving MVP Problem with solve.QP {#exampleanalyticalmvp}
This section provides insights into the effects of diversification and the use of solve.QP by creating ten different efficient frontiers from a pool of ten assets. Each efficient frontier $i \in {1, 2, \cdots, 10}$ consists of $N_i = i$ assets and is created by adding the asset with the next smallest variance first. After loading the returns for ten of the largest stocks in the U.S. market, the variance is calculated to rank all columns in ascending order of variance, as shown in the code below:
returns_raw <- buffer( get_yf( tickers = c("IBM", "GOOG", "AAPL", "MSFT", "AMZN", "NVDA", "JPM", "META", "V", "WMT"), from = "2018-01-01", to = "2019-12-31" )$returns, "AS_10_assets" ) # re-arrange: low var first vars <- sapply(returns_raw, var) returns_raw <- returns_raw[, order(vars, decreasing = F)]
The next step is to create a function mvp() that has the arguments return and lambda. It computes the expected returns mu and the estimated positive definite covariance cov. It then solves an MVP problem with constraints $\textstyle\sum w_i = 1$ and $w_i \geq 0$, which yields the key features mu, var and composition of the portfolio.
mvp <- function(returns, lambda){ tc <- tryCatch({ mu <- ret_to_geomeanret(returns) cov <- as.matrix(nearPD(cov_(returns, mu))$mat) mat <- list( Dmat = lambda * cov, dvec = (1-lambda) * mu, Amat = t(rbind( rep(1, ncol(returns)), # sum up to 1 diag( 1, nrow=ncol(returns), ncol=ncol(returns) ) # long only )), bvec = c( 1, # sum up to 1 rep(0, ncol(returns)) # long only ), meq = 1 ) qp <- solve.QP( Dmat = mat$Dmat, dvec = mat$dvec, Amat = mat$Amat, bvec = mat$bvec, meq = mat$meq ) res <- list( "mu" = mu %*% qp$solution, "var" = t(qp$solution) %*% cov %*% qp$solution, "composition" = setNames(qp$solution, colnames(returns)) ) TRUE }, error = function(e){FALSE}) if(tc){ return(res) }else{ return(list( "mu" = NA, "var" = NA, "composition" = NA )) } }
Each $\lambda \in {0.01, 0.02, \cdots, 1}$ and each combination of ascending number of assets results in a portfolio that can be created with two for loops.
df <- data.frame( "index"=1, "var"=as.numeric(var(returns_raw[, 1])), "return" = as.numeric(ret_to_geomeanret(returns_raw[, 1])), row.names=NULL ) for(i in 2:ncol(returns_raw)){ returns <- returns_raw[, 1:i] for(lambda in seq(0.01, 1, 0.01)){ res <- mvp(returns, lambda) df <- rbind( df, data.frame("index"=i, "var"=res$var, "return" = res$mu) ) } }
The result is filtered and names are added to represent the number of assets. Now the $\mu$-$\sigma$ diagram can be created:
df <- df %>% filter(!is.na(return)) %>% distinct() %>% mutate(name = paste0("n_", index)) %>% arrange(name) %>% mutate(name = factor(name, levels=paste0("n_", ncol(returns_raw):1))) max_show_sd <- df %>% group_by(index) %>% summarise(max_x = max(var)) %>% pull(max_x) %>% mean() %>% sqrt() plot_ly( data = df[df$index!=1,], x=~sqrt(var), y=~return, name=~name, mode="lines", type = 'scatter', color = ~name, colors = c("green", "red") ) %>% add_trace( data=df[df$index==1,], x=~sqrt(var), y=~return, showlegend=T, marker=list(color="red"), mode="markers", name="n_1") %>% layout( xaxis=list(range=c(sqrt(min(df$var))*0.9, max_show_sd), title="sigma"), yaxis=list(range=c(min(df$return)*0.9, (max(df$return)+mean(df$return))*0.5), title="mu")) %>% html_save()
It can be seen that each added asset leads to a minimum variance portfolio with a smaller or equal standard deviation. Despite starting with the asset that has the smallest standard deviation of r sqrt(min(vars)). This is the effect of diversification mentioned by Markowitz.
Solving ITP-MSTE with solve.QP {#exampleitpsolveqp}
This example analyzes how many assets are needed to minimize the mean square error between the replication and historical returns of the SP500TR from 2018-01-01 to 2019-12-31. The constraints are set to be long only and the weights should sum to one. To gradually reduce the number of assets, the five assets with the lowest weights are discarded, and the remaining assets serve as the new asset pool for the next replication until only five assets remain. First, the required data can be downloaded from the R/ directory using existing functions. The function get_spx_composition() uses web scraping to read the components of wikipedia and converts them into monthly compositions of the SP500TR. The pool is formed from all assets present in the last month of the time frame, reduced by assets with missing values. The code below loads the returns of all assets in the pool and the SP500TR:
from <- "2018-01-01" to <- "2019-12-31" spx_composition <- buffer( get_spx_composition(), "AS_spx_composition" ) pool_returns_raw <- buffer( get_yf( tickers = spx_composition %>% filter(Date<=to) %>% filter(Date==max(Date)) %>% pull(Ticker), from = from, to = to )$returns, "AS_sp500_assets" ) pool_returns_raw <- pool_returns_raw[, colSums(is.na(pool_returns_raw))==0] bm_returns <- buffer( get_yf(tickers = "^SP500TR", from = from, to = to)$returns, "AS_sp500tr" ) %>% setNames(., "SP500TR")
The required data is now available and the function for the ITP-MSTE can be created. It requires pool_returns with variable number of columns and the single-column matrix bm_returns.
itp <- function(pool_returns, bm_returns){ mat <- list( Dmat = t(pool_returns) %*% pool_returns, dvec = t(pool_returns) %*% bm_returns, Amat = t(rbind( rep(1, ncol(pool_returns)), # sum up to 1 diag(1, nrow=ncol(pool_returns), ncol=ncol(pool_returns)) # long only )), bvec = c( 1, # sum up to 1 rep(0, ncol(pool_returns)) # long only ), meq = 1 ) qp <- solve.QP( Dmat = mat$Dmat, dvec = mat$dvec, Amat = mat$Amat, bvec = mat$bvec, meq = mat$meq ) res <- list( "var" = as.numeric( var(pool_returns %*% qp$solution - bm_returns)), "solution" = setNames(qp$solution, colnames(pool_returns)) ) }
The successive removal of assets can begin and the results are stored in res.
res <- NULL asset_pool <- NULL n_assets <- rev(seq(5, ncol(pool_returns_raw), 5)) for(i in n_assets){ temp <- if(i==max(n_assets)){ itp(pool_returns_raw, bm_returns) }else{ asset_pool <- names(sort(temp$solution, decreasing = T)[1:i]) itp( pool_returns_raw[, asset_pool], bm_returns ) } res <- rbind( res, data.frame("N"=i, "var"=temp$var, "sd"=sqrt(temp$var), row.names = NULL) ) # for later analysis if(length(asset_pool)==100){ assets_pool_100 <- asset_pool save(assets_pool_100, file="data/assets_pool_100.rdata") } if(length(asset_pool)==50){ assets_pool_50 <- asset_pool save(assets_pool_50, file="data/assets_pool_50.rdata") } }
The following chart illustrates the increase in standard deviation to track the historical performance of the SP500TR as the number of assets $N$ in the asset pool decreases:
plot_ly(data=res, x=~N, y=~sd, mode="lines", type = 'scatter') %>% layout(yaxis=list(range=c(0, mean(max(res$sd),mean(res$sd)) ))) %>% html_save()
It can be seen that the standard deviation stagnates at about $N=200$. This leads to the conclusion that a sparse replication with two hundred assets is sufficient in this particular case to track the historical performance of the SP500TR over this period.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.