R中的递归优化（）函数会导致错误

Question

我正在尝试使用 R 中的 optim() 函数来最小化矩阵运算的值。在这种情况下，我试图最小化一组股票的波动性，这些股票的个人回报彼此共变。被最小化的目标函数是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))

当我手动调用目标函数时，它会按预期返回一个值：

> calculate_portfolio_variance(allocations, returns_covar, ticker_avg, benchmark_avg)
          [,1]
[1,] 0.1713439

但是，当我使用 optim() 函数时，它会返回错误：

> 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

我不确定原因，但可能与optim() 递归使用allocations 变量的方式有关。我能做些什么来解决这个问题？

编辑：FWIW，其他优化策略有效（差分进化，模拟退火），但我更喜欢使用梯度下降，因为它要快得多

Answer 1

如果将第一个参数重命名为 par 并且您切换将 t() 应用于该侧翼矩阵乘法运算中使用的参数向量的顺序，则不会发生错误：

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)
}

我在代码中留下了par的调试打印，并显示了运行结果的顶部

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.

... 结果：

> 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"

正如我在对不相关问题的评论中所说，优化函数首先尝试提高然后降低 par 中的第一个元素，然后尝试对第二个、第三个和第四个元素执行相同的操作。那时发现没有任何改进，它“决定”它收敛到局部最小值并宣布收敛。

我应该指出optim 的代码是rather old and the author of the original algorithm, Dr Nash，已经以the optimx package 的形式在CRAN 上放置了一个更新版本。他说optim 在当时很好，但他认为如果不成功，应该尝试其他程序。