Portfolio Optimization and Returns in R

Implementation in R

       For this task, we will import our data from Yahoo Finance. The five assets under examination are Exchange Traded Funds, which are funds that can be traded on an exchange like a stock. Exchange-traded funds are a type of investment fund that offers the best attributes of two popular assets: they have the diversification benefits of mutual funds and the ease with which stocks are traded. However, before we import any data, we must answer the following question: In what form do we want our data to be? Since we are in the world of finance, times series is the most common type of data. R has the xts package to handle data that are indexed by date. Our task, therefore, reduces to the following:

• Import daily prices from Yahoo Finance

• Convert daily prices to monthly return

• Because we will be aggregating monthly returns to form a portfolio, we will need to compute the simple returns

# Create a vector of ticker symbols
symbols <- c("SPY", "EFA", "IJS", "EEM", "AGG")
# Load data from 2012 to today
# Specify the "to = " argument to specify an end date
prices <- quantmod::getSymbols(
  Symbols = symbols,
  src = "yahoo",
  from = "2012-12-31",
  auto.assign = TRUE,
  warnings = FALSE
) %>%
  # The map function takes an anonymous function and will return a list of five
  # The function Ad() extracts the daily adjusted price series
  purrr::map(.f = ~ quantmod::Ad(get(x = .x))) %>%
  # Use reduce() to merge the elements of .x interactively
  purrr::reduce(.f = merge) %>%
  # Use a replacement function to set column names to ticker symbols
  # This function is in prefix form
  # It is equivalent to colnames(x = prices) <- value
  `colnames<-`(value = symbols)
# Keep only the last reading of each month
# We could have chosen to keep only the first reading of each month
asset_returns_xts <- xts::to.monthly(
  x = prices,
  drop.time = TRUE,
  indexAt = "lastof",
  OHLC = FALSE
) %>%
  # Compute simple returns
  # Log returns are time-additive but not portfolio additive
  PerformanceAnalytics::Return.calculate(method = "discrete") %>%
  # Drop the first row since we lose 12/31/2012
  stats::na.omit()
# Keep only the xts returns, ticker symbols, and the prices series
rm(list = setdiff(x = ls(), y = c("symbols", "prices", "asset_returns_xts")))

• Create Portfolio object

# Examine the monthly simple returns for our five ETF's
head(x = asset_returns_xts, 5)

                  SPY         EFA         IJS          EEM           AGG
2013-01-31 0.05119052  0.03728466 0.053516334 -0.002930946 -0.0062115144
2013-02-28 0.01275821 -0.01288577 0.016306731 -0.022840526  0.0059086990
2013-03-31 0.03797176  0.01305398 0.041079862 -0.010182972  0.0009852088
2013-04-30 0.01921213  0.05018653 0.001223168  0.012158139  0.0096857020
2013-05-31 0.02360972 -0.03019051 0.042869516 -0.048279085 -0.0200111498

# Create Portfolio object which is essentially a list object
min_var_portfolio <- PortfolioAnalytics::portfolio.spec(assets = symbols)
typeof(min_var_portfolio)

[1] "list"

• Add constraints to the portfolio object

# Add the full investment constraint that specifies that the weights must sum to 1
min_var_portfolio <- PortfolioAnalytics::add.constraint(
  portfolio = min_var_portfolio,
  type = "full_investment"
)
# Examine the constraint element by extracting min_var_portfolio[["constraints"]][[1]]
str(pluck(.x = min_var_portfolio, "constraints", 1))

List of 6
 $ type   : chr "full_investment"
 $ enabled: logi TRUE
 $ message: logi FALSE
 $ min_sum: num 1
 $ max_sum: num 1
 $ call   : language PortfolioAnalytics::add.constraint(portfolio = min_var_portfolio, type = "full_investment")
 - attr(*, "class")= chr [1:2] "weight_sum_constraint" "constraint"

# Add the box constraint that ensure the weights are between 0.1 and 0.6
min_var_portfolio <- PortfolioAnalytics::add.constraint(
  portfolio = min_var_portfolio,
  type = "box", min = 0.05, max = 0.6
)
# Examine the constraint element by extracting min_var_portfolio[["constraints"]][[2]]
str(pluck(.x = min_var_portfolio, "constraints", 2))

List of 5
 $ type   : chr "box"
 $ enabled: logi TRUE
 $ min    : Named num [1:5] 0.05 0.05 0.05 0.05 0.05
  ..- attr(*, "names")= chr [1:5] "SPY" "EFA" "IJS" "EEM" ...
 $ max    : Named num [1:5] 0.6 0.6 0.6 0.6 0.6
  ..- attr(*, "names")= chr [1:5] "SPY" "EFA" "IJS" "EEM" ...
 $ call   : language PortfolioAnalytics::add.constraint(portfolio = min_var_portfolio, type = "box",      min = 0.05, max = 0.6)
 - attr(*, "class")= chr [1:2] "box_constraint" "constraint"

• Add objective function

# Add objective to minimize variance
min_var_portfolio <- PortfolioAnalytics::add.objective(
  portfolio = min_var_portfolio,
  # Minimize risk
  type = "risk",
  # A character corresponding to a function name, var()
  name = "var"
)

• Optimization

# Run the optimization
global_min_portfolio <- PortfolioAnalytics::optimize.portfolio(
  R = asset_returns_xts,
  portfolio = min_var_portfolio,
  # This defaults to the "quadprog" solver
  optimize_method = "quadprog",
  # Return additional information on the path or portfolios searched
  trace = TRUE
)
# Examine returned portfolio list object
global_min_portfolio

***********************************
PortfolioAnalytics Optimization
***********************************

Call:
PortfolioAnalytics::optimize.portfolio(R = asset_returns_xts, 
    portfolio = min_var_portfolio, optimize_method = "quadprog", 
    trace = TRUE)

Optimal Weights:
  SPY   EFA   IJS   EEM   AGG 
0.213 0.087 0.050 0.050 0.600 

Objective Measure:
 StdDev 
0.01644 

Implementation in R

• Create portfolio object

# Create Portfolio object
max_exp_return_portfolio <- PortfolioAnalytics::portfolio.spec(assets = symbols)

• Add constraints to the object

# Add the full investment constraint that specifies the weights must sum to 1
max_exp_return_portfolio <- PortfolioAnalytics::add.constraint(
  portfolio = max_exp_return_portfolio,
  type = "full_investment"
)
# Add the box constraint that ensure the weights are between 0.1 and 0.6
max_exp_return_portfolio <- PortfolioAnalytics::add.constraint(
  portfolio = max_exp_return_portfolio,
  type = "box", min = 0.05, max = 0.6
)

• Add objective function

# Add objective to maximize mean returns
max_exp_return_portfolio <- PortfolioAnalytics::add.objective(
  portfolio = max_exp_return_portfolio,
  # Maximize expected returns
  type = "return",
  # A character corresponding to a function name, mean()
  name = "mean"
)

• Optimization

# Run the optimization
global_max_portfolio <- PortfolioAnalytics::optimize.portfolio(
  R = asset_returns_xts,
  portfolio = max_exp_return_portfolio,
  # This defaults to the "glpk" solver
  optimize_method = "glpk",
  # Return additional information on the path or portfolios searched
  trace = TRUE
)
# Examine returned portfolio list object
global_max_portfolio

***********************************
PortfolioAnalytics Optimization
***********************************

Call:
PortfolioAnalytics::optimize.portfolio(R = asset_returns_xts, 
    portfolio = max_exp_return_portfolio, optimize_method = "glpk", 
    trace = TRUE)

Optimal Weights:
 SPY  EFA  IJS  EEM  AGG 
0.60 0.05 0.25 0.05 0.05 

Objective Measure:
  mean 
0.0114 

Monthly portfolio returns by hand

       We can compute portfolio monthly returns using a brute force method:

# Compute by hand
portfolio_returns_by_hand <-
  (weights[[1]] * asset_returns_xts[, 1]) +
  (weights[[2]] * asset_returns_xts[, 2]) +
  (weights[[3]] * asset_returns_xts[, 3]) +
  (weights[[4]] * asset_returns_xts[, 4]) +
  (weights[[5]] * asset_returns_xts[, 5])
# Name the series
portfolio_returns_by_hand <- `names<-`(portfolio_returns_by_hand, "Monthly portfolio returns")
# Examine
head(portfolio_returns_by_hand, 5)

           Monthly portfolio returns
2013-01-31                0.04550050
2013-02-28                0.01024073
2013-03-31                0.03324583
2013-04-30                0.01543459
2013-05-31                0.01995917

Monthly portfolio returns in xts

       Another way to compute portfolio monthly returns is to use functions from the PerformanceAnalytics package, which works well with xts objects. We also adopt monthly re-balancing as a strategy. The re-balancing of investments is the action of keeping the portfolio weights consistent with the optimal weights. Note that re-balancing the weights every month may be unrealistic; in reality, the convention is often to re-balance quarterly or annually. For this example, however, we will re-balance monthly to be consistent with our brute force computation.

# Compute monthly portfolio returns
portfolio_returns_xts_rebalanced_monthly <-
  PerformanceAnalytics::Return.portfolio(
    R = asset_returns_xts,
    weights = weights,
    # Monthly re-balancing
    reblance_on = "months",
    # Use simple/arithmetic chaining to aggregate returns
    geometric = FALSE
  ) %>%
  `colnames<-`("Monthly_portfolio_returns")
# Examine
head(portfolio_returns_xts_rebalanced_monthly, 5)

           Monthly_portfolio_returns
2013-01-31                0.04550050
2013-02-28                0.01024073
2013-03-31                0.03324583
2013-04-30                0.01543459
2013-05-31                0.01995917

• The function Return.portfolio(R, weights = NULL, wealth.index = FALSE, contribution = FALSE, geometric = TRUE, rebalance_on = c(NA, "years", "quarters", "months", "weeks", "days"), value = 1, verbose = FALSE, ...) calculates the returns of a portfolio with the same periodicity of the returns data.

Monthly portfolio returns in the tidyverse

       The tidyverse is a collection of R packages designed with the same underlying philosophy, grammar, and data structures. Simply put, the “tidy” data structure that works well with tidyverse functions is one where every row is an observation and every column is a variable. If we re-examine the xts object called “asset_returns_xts,” we see that every column is a returns series for a particular asset. This is inconsistent with the tidyverse data structure, and so we must convert the xts object to a tidy format for our computations. Ideally, we would like to have one column called “asset” that specifies the names of the ETF instead of having five individual columns of returns data. This idea will become clearer once we convert our xts object to a tibble.

Now, examine the “tidy” data structure and compare it to the xts object created earlier:

FSA::headtail(asset_returns_long, n = 5)

          date asset       returns
1   2013-01-31   SPY  0.0511905160
2   2013-01-31   EFA  0.0372846569
3   2013-01-31   IJS  0.0535163338
4   2013-01-31   EEM -0.0029309465
5   2013-01-31   AGG -0.0062115144
521 2021-09-30   SPY -0.0224997875
522 2021-09-30   EFA -0.0017360987
523 2021-09-30   IJS -0.0316548170
524 2021-09-30   EEM -0.0198435604
525 2021-09-30   AGG -0.0001984158

Since our returns data is in a “tidy” format, computing portfolio monthly returns is very a straight forward. For those who are familiar with SQL, we are essentially using the CASE statement here.

# Use vectorized if else statements to assign weights according to the asset column
potfolio_returns_dplyr_byhand <- asset_returns_long %>%
  group_by(asset) %>%
  mutate(
    weights = dplyr::case_when(
      asset == symbols[[1]] ~ weights[[1]],
      asset == symbols[[2]] ~ weights[[2]],
      asset == symbols[[3]] ~ weights[[3]],
      asset == symbols[[4]] ~ weights[[4]],
      asset == symbols[[5]] ~ weights[[5]]
    ),
    weighted_returns = weights * returns
  ) %>%
  # Group by date
  # We need to group by date so that the aggregate sum() function is carried out by month
  # For each date, the original series has 5 weighted returns, one for each ETF
  # The results should be 1 portfolio return (the sum of the 5 weighted returns) for each month
  dplyr::group_by(date) %>%
  # Compute monthly portfolio returns
  dplyr::summarize(Monthly_portfolio_returns = sum(weighted_returns))
# Examine the data
head(potfolio_returns_dplyr_byhand, 5)

# A tibble: 5 × 2
  date       Monthly_portfolio_returns
  <date>                         <dbl>
1 2013-01-31                    0.0455
2 2013-02-28                    0.0102
3 2013-03-31                    0.0332
4 2013-04-30                    0.0154
5 2013-05-31                    0.0200

• The function summarize(.data, ..., .groups = NULL) creates a new data frame. It will have one (or more) rows for each combination of grouping variables.

Monthly portfolio returns in tidyquant

       The tidyquant package gives us two ways to compute the portfolio monthly returns.

# Use function from the tidyquant
portfolio_returns_tq_rebalanced_monthly_method_1 <- asset_returns_long %>%
  tidyquant::tq_portfolio(
    assets_col = asset,
    returns_col = returns,
    weights = weights,
    col_rename = "Monthly_portfolio_returns",
    rebalance_on = "months"
  )
# Examine
head(portfolio_returns_tq_rebalanced_monthly_method_1, 5)

# A tibble: 5 × 2
  date       Monthly_portfolio_returns
  <date>                         <dbl>
1 2013-01-31                    0.0455
2 2013-02-28                    0.0102
3 2013-03-31                    0.0332
4 2013-04-30                    0.0154
5 2013-05-31                    0.0200

• The function tq_portfolio(data, assets_col, returns_col, weights = NULL, col_rename = NULL, ...) aggregates a group of returns by asset into portfolio returns.

       An alternative way to load data is to use the tidyquant wrapper function, which automatically coerces the returns into a long tidy format:

# Load data
asset_returns_tq <- tidyquant::tq_get(
  x = symbols,
  get = "stock.prices",
  from = "2012-12-31"
) %>%
  # The asset column is named symbol by default (see body(tidyquant::tq_get))
  dplyr::group_by(symbol) %>%
  # Select adjusted daily prices
  tidyquant::tq_transmute(
    select = adjusted,
    col_rename = "returns",
    # This function is from quantmod
    mutate_fun = periodReturn,
    # These arguments are passed along to the mutate function quantmod::periodReturn
    period = "monthly",
    # Simple returns
    type = "arithmetic",
    # Do not return leading period data
    leading = FALSE,
    # This argument is passed along to xts::to.period, which is wrapped by quantmod::periodReturn
    # We use the last reading of each month to find percentage changes
    indexAt = "lastof"
  ) %>%
  dplyr::rename(asset = symbol) %>%
  na.omit()

• The function tq_get(x, get = "stock.prices", complete_cases = TRUE, ...) gets data in tibble format. The most important argument is perhaps the dot-dot-dot, since these are the arguments passed to the underlying functions from the other packages that tq_get() uses. There could be multiple layers of wrapper functions, and so it is best practice to read the documentations carefully.

       A possibly more useful method of aggregating asset returns to portfolio returns is then to map a tibble of ticker symbols and weights to the tibble of ticker symbols and monthly returns:

# Create a tibble of weights
weights_tibble <- tibble(
  asset = symbols,
  weights = weights
)
head(weights_tibble, 5)

# A tibble: 5 × 2
  asset weights
  <chr>   <dbl>
1 SPY     0.6  
2 EFA     0.05 
3 IJS     0.250
4 EEM     0.05 
5 AGG     0.05 

# Map the weights to the returns column using the asset column as the identifier
portfolio_returns_tq_rebalanced_monthly_method_2 <- asset_returns_tq %>%
  tidyquant::tq_portfolio(
    assets_col = asset,
    returns_col = returns,
    weights = weights_tibble,
    col_rename = "Monthly_portfolio_returns",
    rebalance_on = "months"
  )
head(portfolio_returns_tq_rebalanced_monthly_method_2, 5)

# A tibble: 5 × 2
  date       Monthly_portfolio_returns
  <date>                         <dbl>
1 2013-01-31                    0.0455
2 2013-02-28                    0.0102
3 2013-03-31                    0.0332
4 2013-04-30                    0.0154
5 2013-05-31                    0.0200

Compare and Contrast the four methods

       We have covered a variety of different methods for aggregating asset monthly returns to portfolio monthly returns. As a sanity check, we want to ensure that these methods yield consistent results:

potfolio_returns_dplyr_byhand %>%
  # Rename the column
  rename(tidyverse_method = Monthly_portfolio_returns) %>%
  # Create three new columns that contain results for other methods
  dplyr::mutate(
    "equation_byhand" = zoo::coredata(x = portfolio_returns_by_hand),
    "tq_method_1" = portfolio_returns_tq_rebalanced_monthly_method_1[["Monthly_portfolio_returns"]],
    "tq_method_2" = portfolio_returns_tq_rebalanced_monthly_method_2[["Monthly_portfolio_returns"]],
    "xts_method" = zoo::coredata(portfolio_returns_xts_rebalanced_monthly)
  ) %>%
  purrr::modify_if(.p = is.numeric, .f = ~ round(x = .x, digits = 7)) %>%
  dplyr::select(-date) %>%
  # Examine
  head(10)

# A tibble: 10 × 5
   tidyverse_method equation_byhand[,… tq_method_1 tq_method_2 xts_method[,"Mon…
              <dbl>              <dbl>       <dbl>       <dbl>             <dbl>
 1           0.0455             0.0455      0.0455      0.0455            0.0455
 2           0.0102             0.0102      0.0102      0.0102            0.0102
 3           0.0332             0.0332      0.0332      0.0332            0.0332
 4           0.0154             0.0154      0.0154      0.0154            0.0154
 5           0.0200             0.0200      0.0200      0.0200            0.0200
 6          -0.0131            -0.0131     -0.0131     -0.0131           -0.0131
 7           0.0509             0.0509      0.0509      0.0509            0.0509
 8          -0.0292            -0.0292     -0.0292     -0.0292           -0.0292
 9           0.0436             0.0436      0.0436      0.0436            0.0436
10           0.0406             0.0406      0.0406      0.0406            0.0406

• The function coredata(x, ...) is a generic function for extracting the core data contained in, a (more complex) object and replacing it. The replacement function is coredata<-.

       As can be seen, the results are all consistent. Sigh of relief. For the upcoming posts, I also plan on tackling other topics in portfolio theory and to build off of what I’ve covered in this post. In particular, Python also has many great libraries— NumPy, SciPy, Quandl— all of which contain great tools for financial analysis. Utilizing R’s tidyverse ecosystem and Python’s fast, memory-efficient methods, there’s great deal of content to cover.