I am trying to use the optim() function in R to minimize a value with matrix operations. In this case, I am trying to minimize the volatility of a group of stocks whose individual returns covary with each other. The objective function being minimized is calculate_portfolio_variance.
library(quantmod)
filter_and_sort_symbols <- function(symbols)
{
# Name: filter_and_sort_symbols
# Purpose: Convert to uppercase if not
# and remove any non valid symbols
# Input: symbols = vector of stock tickers
# Output: filtered_symbols = filtered symbols
# convert symbols to uppercase
symbols <- toupper(symbols)
# Validate the symbol names
valid <- regexpr("^[A-Z]{2,4}$", symbols)
# Return only the valid ones
return(sort(symbols[valid == 1]))
}
# Create the list of stock tickers and check that they are valid symbols
tickers <- filter_and_sort_symbols(c("AAPL", "NVDA", "MLM", "AA"))
benchmark <- "SPY"
# Set the start and end dates
start_date <- "2007-01-01"
end_date <- "2019-01-01"
# Gather the stock data using quantmod library
getSymbols(Symbols=tickers, from=start_date, to=end_date, auto.assign = TRUE)
getSymbols(benchmark, from=start_date, to=end_date, auto.assign = TRUE)
# Create a matrix of only the adj. prices
price_matrix <- NULL
for(ticker in tickers){price_matrix <- cbind(price_matrix, get(ticker)[,6])}
# Set the column names for the price matrix
colnames(price_matrix) <- tickers
benchmark_price_matrix <- NULL
benchmark_price_matrix <- cbind(benchmark_price_matrix, get(benchmark)[,6])
# Compute log returns
returns_matrix <- NULL
for(ticker in tickers){returns_matrix <- cbind(returns_matrix, annualReturn(get(ticker)))}
returns_covar <- cov(returns_matrix)
colnames(returns_covar) <- tickers
rownames(returns_covar) <- tickers
# get average returns for tickers and benchmark
ticker_avg <- NULL
for(ticker in tickers){ticker_avg <- cbind(ticker_avg, colMeans(annualReturn(get(ticker))))}
colnames(ticker_avg) <- tickers
benchmark_avg <- colMeans(annualReturn(get(benchmark)))
# create the objective function
calculate_portfolio_variance <- function(allocations, returns_covar, ticker_avg, benchmark_avg)
{
# Name: calculate_portfolio_variance
# Purpose: Computes expected portfolio variance, to be used as the minimization objective function
# Input: allocations = vector of allocations to be adjusted for optimality; returns_covar = covariance matrix of stock returns
# ticker_avg = vector of average returns for all tickers, benchmark_avg = benchmark avg. return
# Output: Expected portfolio variance
# get benchmark volatility
benchmark_variance <- (sd(annualReturn(get(benchmark))))^2
# scale allocations for 100% investment
allocations <- as.matrix(allocations/sum(allocations))
# get the naive allocations
naive_allocations <- rep(c(1/ncol(ticker_avg)), times=ncol(ticker_avg))
portfolio_return <- sum(t(allocations)*ticker_avg)
portfolio_variance <- t(allocations)%*%returns_covar%*%allocations
# constraints = portfolio expected return must be greater than benchmark avg. return and
# portfolio variance must be less than benchmark variance (i.e. a better reward at less risk)
if(portfolio_return < benchmark_avg | portfolio_variance > benchmark_variance)
{
allocations <- naive_allocations
}
portfolio_variance <- t(allocations)%*%returns_covar%*%allocations
return(portfolio_variance)
}
# Specify lower and upper bounds for the allocation percentages
lower <- rep(0, ncol(returns_matrix))
upper <- rep(1, ncol(returns_matrix))
# Initialize the allocations by evenly distributing among all tickers
set.seed(1234)
allocations <- rep(1/length(tickers), times=length(tickers))
When I call the objective function manually, it returns a value as expected:
> calculate_portfolio_variance(allocations, returns_covar, ticker_avg, benchmark_avg)
[,1]
[1,] 0.1713439
However, when I use the optim() function it returns the error:
> optim_result <- optim(par=allocations, fn=calculate_portfolio_variance(allocations, ticker_avg, benchmark_avg), lower=lower, upper=upper, method="L-BFGS-B")
Error in t(allocations) %*% returns_covar : non-conformable arguments
I'm not sure the reason, but it may be with the way optim() recursively uses the allocations variable. What can I do to fix this?
Edit: FWIW, other optimization strategies work (differential evolution, simulated annealing) but I would prefer to use gradient descent because it is considerably faster
No error occurs if the first argument is renamed par and you switch the order in which you apply t() to the parameter vectors used in that flanking matrix-multiply operation:
cpv <- function(par, returns_covar=returns_covar, ticker_avg, benchmark_avg)
{
# Name: calculate_portfolio_variance
# Purpose: Computes expected portfolio variance, to be used as the minimization objective function
# Input: allocations = vector of allocations to be adjusted for optimality; returns_covar = covariance matrix of stock returns
# ticker_avg = vector of average returns for all tickers, benchmark_avg = benchmark avg. return
# Output: Expected portfolio variance
# get benchmark volatility
benchmark_variance <- (sd(annualReturn(get(benchmark))))^2
# scale allocations for 100% investment
par <- as.matrix(par/sum(par))
# get the naive allocations
naive_allocations <- rep(c(1/ncol(ticker_avg)), times=ncol(ticker_avg))
portfolio_return <- sum(t(par)*ticker_avg);print(par)
portfolio_variance <- t(par)%*%returns_covar%*%par
# constraints = portfolio expected return must be greater than benchmark avg. return and
# portfolio variance must be less than benchmark variance (i.e. a better reward at less risk)
if(portfolio_return < benchmark_avg | portfolio_variance > benchmark_variance)
{
par <- naive_allocations
}
portfolio_variance <- t(par)%*%returns_covar%*%par
return(portfolio_variance)
}
I left the debugging printing of par in the code and show the top of the results of running it
optim_result <- optim(par=allocations, fn=cpv, lower=lower, upper=upper, returns_covar=returns_covar, ticker_avg=ticker_avg, benchmark_avg=benchmark_avg, method="L-BFGS-B")
[,1]
[1,] 0.25
[2,] 0.25
[3,] 0.25
[4,] 0.25
[,1]
[1,] 0.2507493
[2,] 0.2497502
[3,] 0.2497502
[4,] 0.2497502
[,1]
[1,] 0.2492492
[2,] 0.2502503
[3,] 0.2502503
[4,] 0.2502503
#--- snipped output of six more iterations.
... and the result:
> optim_result
$par
[1] 0.25 0.25 0.25 0.25
$value
[1] 0.1713439
$counts
function gradient
1 1
$convergence
[1] 0
$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
As I said in the comment to an unrelated question, the optim function first tried to raise then lower the first element in par, then tries to do the same for the second, third and fourth elements. At that point finding no improvement it "decides" it's converged at a local minimum and declares convergence.
I should point out that the code for optim is rather old and the author of the original algorithm, Dr Nash, has placed an updated version on CRAN in the form of the optimx package. He says optim was good in its day, but that he thinks other procedures should be tried if it's not successful.