代码之家  ›  专栏  ›  技术社区  ›  coolhand

R中的递归optim()函数导致错误

  •  0
  • coolhand  · 技术社区  · 5 年前

    我正在尝试使用 optim() R中的函数,通过矩阵运算使值最小化。在这种情况下,我试图最小化一组股票的波动性,这些股票的个人回报是相互协变的。被最小化的目标函数是 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,其他优化策略也有效(差分进化、模拟退火),但我更喜欢使用梯度下降法,因为它要快得多

    0 回复  |  直到 5 年前
        1
  •  1
  •   IRTFM    5 年前

    如果将第一个参数重命名为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"
    

    正如我在对一个无关问题的评论中所说,optim函数首先尝试提高然后降低par中的第一个元素,然后尝试对第二、第三和第四个元素做同样的事情。在这一点上没有发现任何改进,它“决定”它收敛到局部最小值并宣布收敛。

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