This is the supplement to Chapter 4 of my book, Goals-Based Portfolio Theory, demonstrating the techniques and offering some code examples. If you are reading this having not purchased/read the book, you are missing much of the narrative. You can pick up a copy from Wiley, Amazon.com, or Barnes & Noble.

Most of the time, as practitioners, we have a multi-period view of markets and investing. That is, rather than generate long-term capital market assumptions (e.g. SPX yielding an 8% return with 16% volatility over the coming decade), we may wish to express a view that a recession is forming this year, with a recovery the year after, then a return to the long-run average the year after that.

Unfortunately, the standard optimizers in the industry do not allow for a multi-period view. Here we will explore a simple example of how to optimize a goals-based portfolio across multiple periods (for a fuller treatment of the math and theory, you’ll have to buy my book!).

Here is a basic outline of the approach:

1. Develop capital market expectations for each period (in our example, we will consider three periods into the future).
2. Sample each period’s portfolio distribution (which will be determined by the weights of investments in that period) and generate potential “paths” for the portfolio to follow through time.
3. Combine these “paths” to generate an estimate of the probability of failing to achieve the annual return requirement required by the goal.

First, let us load our packages and build the necessary functions.

library(tidyverse)
library(nloptr)
# Kronecker product of a list of vectors, relies on (A %x% B) %x% C = A %x% (B %x% C)
list_kronecker.f <- function( list_of_vectors ){
  proxy <- list( list_of_vectors[[1]] )
  for(i in 2:length(list_of_vectors)){
    proxy[[i]] <- proxy[[i-1]] %x% list_of_vectors[[i]]
  }
  return( proxy[[ length(proxy) ]] )
}

# Expected volatility function
vol.f <- function(weight_vector, covariance){
  sqrt( t(weight_vector) %*% covariance %*% weight_vector )
}

prob.f <- function(x, m, s, conversion_factor){
  1/(1 + exp(-(x-m)/(conversion_factor * s)))
}

# Expected lower-tail cumulative distribution function
list_cdf.f <- function( required_return, # annual return required by the investor
                        returns_list, # asset return expectations by period
                        weights_list, # weights vector list
                        covariance_table_list, # asset covariance expectations by period
                        point_estimates = 60 ){ # number of point estimates
  
  num_assets <- length( returns_list[[1]] ) # number of assets
  num_periods <- length( returns_list ) # number of periods
  xmin <- -0.80 # minimum point estimate in return distribution
  xmax <- 0.80 # maximum point estimate in return distribution
  dx <- (xmax - xmin)/point_estimates # resolution of estimate
  x <- seq(xmin, xmax, dx) # return vector
  
  # Take point estimates of each period's distribution
  D <- vector( mode = 'list', length = num_periods ) # return distribution point estimates
  P <- vector( mode = 'list', length = num_periods ) # probability point estimates
  for(j in 1:num_periods){
    for(i in 1:point_estimates){
      
      D[[j]][i] <- (x[i+1] + x[i])/2 + 1 # return point estimate, add one for later
      m <- sum( weights_list[[j]] * returns_list[[j]] ) # return expectation of these weights
      s <- vol.f( weights_list[[j]], covariance_table_list[[j]] ) # st dev expectation of these weights
      P[[j]][i] <- pnorm(x[i+1], m, s ) - pnorm( x[i], m, s ) # probability point estimate
      # P[[j]][i] <- prob.f(x[i+1], m, s, 0.60) - prob.f(x[i], m, s, 0.60) # For logistic
      
    }
  }
  
  DF <- list_kronecker.f( D )^(1/num_periods) - 1 # Apply geometric return pattern for final distribution
  PF <- list_kronecker.f( P ) # Final probability distribution
  o <- order( DF, decreasing = F ) # Order the returns from low to high
  DF <- DF[ o ]
  PF <- PF[ o ] # align the probability distribution with returns order
  index <- which( DF <= required_return ) # index those returns <= threshold return
  return(
    cumsum(PF)[ index ][ length(index) ] # return the cumsum up to the point of the threshold return
  )
}

# Optimizer function--input Multi-Period CMEs and output optimal weights
MPGB_optimize.f <- function( required_return, # Required annualized return
                             returns_list, # A list of return vectors
                             covariance_table_list, # A list of covariance matrices
                             point_estimates = 60 ){ # number of point estimates
  
  num_assets <- length( returns_list[[1]] ) # number of assets
  num_periods <- length( returns_list ) # number of periods
  
  # Build constraint function, ensure weights in each period sum to 1
  constraint_function.f <- function( weights_vector ){
    
    const <- 0
    for(i in 1:num_periods){
      end <- num_assets * i
      start <- end - num_assets + 1
      const[i] <- (sum( weights_vector[ start:end ] ) - 1)^2
    }
    return(
      sum( const )
    )
  }
  
  # Build optimization function
  optim_function.f <- function( weights_vector ){
    
    # Convert vectorized weights into a list
    weights_list <- vector( mode = 'list', length = num_periods ) 
    for(i in 1:num_periods){
      end <- num_assets * i
      start <- end - num_assets + 1
      weights_list[[i]] <- weights_vector[ start:end ]
    }
    return(
      list_cdf.f( required_return,
                  returns_list,
                  weights_list,
                  covariance_table_list,
                  point_estimates = 60 ) +
        100 * constraint_function.f( weights_vector )
    )
  }

Note that these functions expect our capital market expectations to be in lists. Since CMEs are our first major step, let’s lay out our assets and our expectations for them over the coming periods.

• Asset A: a large-cap stock-like asset.
• Asset B: a gold-like asset
• Asset C: a bond-like asset

And, just to illustrate the point, we are going to assume that our view of the coming periods is like this:

• Period 1: A recession year
• Period 2: A recovery year
• Period 3: An average year

This view is reflected in the listed CMEs here:

# Set capital market expectations at t1, t2, and t3
# Asset 1: SPX-like asset
# Asset 2: Gold-like asset
# Asset 3: Bond-like asset
# t1 represents a "recession" type year, with mean = -12% and sd = 30%
# t2 represents a "recovery" year, with mean = 16%, and  sd = 15%
# t3 represents an "average" year, with mean = 10%, and  sd = 18%

returns_list <- list( c( -0.2, 0.04, 0.03 ),
                      c( 0.12, 0.04, 0.03 ),
                      c( 0.10, 0.05, 0.04 ) )

# Covar table is formatted as correlation * sd of column * sd of row
covariance_table_list <- list( 
  matrix( c(1.00*0.30*0.30, -0.10*0.30*0.20, 0.10*0.30*0.05,
            -0.10*0.30*0.20, 1.00*0.15*0.20, 0.00*0.20*0.05,
            0.10*0.30*0.05, 0.00*0.05*0.20, 1.00*0.05*0.05 ),
          nrow = 3, ncol = 3, byrow = T),
  
  matrix( c(1.00*0.15*0.15, 0.10*0.20*0.15, 0.10*0.15*0.03,
            0.10*0.20*0.15, 1.00*0.20*0.20, 0.00*0.20*0.03,
            0.10*0.15*0.03, 0.00*0.20*0.03, 1.00*0.03*0.03),
          nrow = 3, ncol = 3, byrow = T),
  
  matrix( c(1.00*0.18*0.18, 0.10*0.18*0.20, 0.10*0.20*0.03,
            0.10*0.18*0.20, 1.00*0.20*0.20, 0.00*0.20*0.03,
            0.10*0.20*0.03, 0.00*0.20*0.03, 1.00*0.03*0.03),
          nrow = 3, ncol = 3, byrow = T)
  )

And, from here, we find the optimal weights to each investment in each period.

(weights_optimum <- MPGB_optimize.f( 0.07,
                                    returns_list,
                                    covariance_table_list,
                                    point_estimates = 60 ))

Yielding an intuitive optimal allocation across each period. Notice the reliance on gold and bonds through the recession, with an all-in on stocks during the recovery and return to average.

• Period 1 (recessionary)
  • 0% to Asset A,
  • 25% to Asset B, and
  • 75% to Asset C
• Period 2 (recovery)
  • 100% to Asset A,
  • 0% to Asset B, and
  • 0% to Asset C
• Period 3 (average)
  • 100% to Asset A,
  • 0% to Asset B, and
  • 0% to Asset C

This multi-period plan of investing generates a 45% probability of failing to achieve a 7% return on average, over the three years.

We can compare this multi-period approach to the single-period optimization approach (note that we have taken the long-term view of Period 3 in our CMEs for this optimization method).

SP_optim_fun <- function( weights ){
  pnorm(required_return,
        mean = sum( returns_list[[3]] * weights),
        sd = vol.f( weights, covariance_table_list[[3]]) ) +
    100 * (sum(weights) - 1)^2
}
required_return <- 0.07
wx <- runif(3) # Generate weight seeds
w <- wx/sum(wx) # Generate random starting weights

(weights_optimum_single_period <- slsqp( w,
                                        SP_optim_fun,
                                        lower = rep(0, length(w) ),
                                        upper = rep(1, length(w) ),
                                        control = list( ftol_rel = 1e-4,
                                                        xtol_rel = 1e-4) )$par )

Which yields a non-intuitive constant allocation across all periods of

• 100% to Asset A
• 0% to Asset B
• 0% to Asset C

And by applying our custom list_cdf.f() function, we find that this allocation yields a 71% probability of failing to achieve our 7% required return! Quite a difference from the 45% of the multi-period optimizer.

list_cdf.f( 0.07,
            returns_list,
            list( weights_optimum_single_period,
                  weights_optimum_single_period,
                  weights_optimum_single_period ),
            covariance_table_list,
            point_estimates = 60 )

The computational burden of this approach is a headwind, I must admit.[1] Though, a more skilled architect might build a more efficient execuation algorithm using this method. That would be helpful since I more-or-less apply this same method to another common (and largely unsolved) portfolio problem: the combination of non-Gaussian asset return distributions.

But that is a post for another day.

[1] I want to also direct interested readers to other resources on the topic, especially the methods of Das, et al (2022), and Das, et al (2019), both of which are more computationally efficient than the method I propose, a consideration which may be important depending on the application.