*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:

- Develop capital market expectations for each period (in our example, we will consider three periods into the future).
- 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.
- 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.

Franklin,

Thank you for this blog.

One question: the function “MPGB_optimize.f” doesn’t have a closing “}” on the first block of code. What would be the correct place?

Tks

LikeLiked by 1 person

The closing bracket for the function is under the return( sum( const ) ) commands:

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 )

)

}

Is that the one you were looking for? Or is there another that I'm missing?

LikeLike