Portfolio basics

This is the companion notebook to the introductory course on simple portfolio construction and performance metrics. It contains simple examples of portfolio backtesting. The structures of this section will be used later on in the course.

\(\rightarrow\) ABOUT THE NOTEBOOK:

• the code chunks are sequential. They must be executed in the correct order.
  
• text between stars appears in bold in html format.

1 A glimpse at the S&P500

The S&P500 is one of the most widely scrutinised index in the US Equity investment space. It serves very often as benchmark. The purpose of this section is to unveil a few of its statistical properties and to use the highcharter package for highchart rendering in R.

1.1 Prices

We start by loading the packages and downloading the data. For technical reasons, we download the series of the ETF replicating the S&P500 (SPY ticker).

if(!require(highcharter)){install.packages("highcharter")} # Package for nice financial graphs
if(!require(scales)){install.packages("scales")}           # Package for graph scales
library(tidyverse)   # The data wrangling package
library(plotly)      # For interactive graphs
library(lubridate)   # For data management
library(quantmod)    # The package that eases the downloading of financial data
library(highcharter) # The package for financial time-series plots

min_date <- "1980-01-01"
max_date <- "2025-09-30"
prices <- getSymbols("SPY", src = 'yahoo',  # Yahoo source (ETF ticker) ^GSPC
                     from = min_date, 
                     to = max_date,
                     auto.assign = TRUE, 
                     warnings = FALSE) %>%
  map(~Ad(get(.))) %>% 
  reduce(merge) %>%
  `colnames<-`("SPY")

Then, we turn to plotting, using highchart format. We underline that this package works with special xts (R extensible time-series) format. We point to the package reference for more details on this subject:
- https://jkunst.com/highcharter/
- https://www.highcharts.com/blog/tutorials/highcharts-for-r-users/
- https://cran.r-project.org/web/packages/highcharter/index.html

highchart(type = "stock") %>%
  hc_title(text = "Evolution of the SPY") %>%
  hc_add_series(prices) %>%
  hc_tooltip(pointFormat = '{point.y:.3f}') # To round the y-axis numbers

A nice feature of highcharts is that they allow the user to change the observation period and see the values of points on the curve.

1.2 Returns

Next, we turn to the distribution of returns.

returns <- prices %>%                            # Formula for returns
  data.frame(Date = index(.)) %>%              # Adding the Date as a column
  mutate(SPY = SPY / dplyr::lag(SPY) - 1) %>%  # Returns
  na.omit()                                    # Removing NA rows
m <- mean(returns$SPY)                           # Average daily return
s <- sd(returns$SPY)                             # Volatility = sd of daily returns

returns %>% ggplot() +                         # Plot
  geom_histogram(aes(x = SPY, y = ..density..), bins = 100, alpha = 0.6) + theme_light() +
  stat_function(fun = dnorm, args = list(mean = m, sd = s), aes(color = "Gaussian")) +
  theme(legend.position = c(0.7, 0.7))

#ggsave("Hist_SPX.pdf", width = 8, height = 5)

We plot the Gaussian density with parameters corresponding to the sample mean and standard deviation. Small grey rectangles around \(\pm 0.05\) indicate that large positive and negative returns occur more often than estimate by the Gaussian law: the tails of their distribution are notoriously heavy.

dens2 <- function(x) dnorm(x, mean = m, sd = s)/mean(returns<(-0.02) & returns>(-0.12), na.rm=T) # Need to rescale
returns %>% ggplot() +                         # Plot
  geom_histogram(aes(x = SPY, y = ..density..), bins = 100, alpha = 0.6) + theme_light() +
  stat_function(fun = dens2, aes(color = "Gaussian")) + 
  theme(legend.position = c(0.4, 0.4)) + 
  xlim(-0.12,-0.02)

#ggsave("Tail_SPX.pdf", width = 8, height = 5)

The Gaussian density is above the histogram for returns above -3%, but then strongly below and predicts that events such as -5% on a given day should virtually never happen.

1.3 Volatility

Finally, we take a dynamic look at the volatility. We take a frugal approach; much more elegant methods are presented in Section 4.5 of Reproducible Finance by Jonathan Regenstein, as well as the roll_sd() function used in http://www.reproduciblefinance.com/drafts/post/vix-and-realized-volatility-updating-our-previous-work/

nb_days <- 21                                       # 63 days roughly equivalent to 3 months 
vol <- 0                                            # Initialisation
for(i in 1:(nrow(returns) - nb_days + 1)){          # Loop on dates: not elegant!
  vol[i] <- sd(returns$SPY[i:(i + nb_days - 1)])    # Forward vol computed on rolling window of nb_days
}
Date <- returns$Date[nb_days:nrow(returns)]         # Vector of dates
vol <- data.frame(vol*sqrt(252))                    # Vector inside a dataframe

rownames(vol) <- Date                               # Change the index of the dataframe
highchart(type = "stock") %>%                       # Code for the plot
  hc_title(text = "Evolution of realized volatility") %>%
  hc_add_series(as.xts(vol)) %>%
  hc_tooltip(pointFormat = '{point.y:.2f}') # To round the y-axis numbers

Clearly the graph shows periods of low volatility and clusters of high market turbulence (crashes, most of the time). Returns are therefore not stationary. The properties we exhibited for the S&P500 are also true at the individual stock level.

In the next sections, we turn to the core topic of portfolio backtesting.

2 Data

Portfolio backtesting usually involves several assets. Hence we rely on an external dataset henceforth.
First, let’s start with the preprocessing of the data.

We first import and arrange the data. The data consists of monhtly financial information pertaining to 30 large US firms. They are characterised by their ticker symbol:

A - F	G - M	O - Z
AAPL (Apple)	GE (General Electric)	ORCL (Oracle)
BA (Boeing)	HD (Home Depot)	PFE (Pfizer)
BAC (Bank of America)	IBM	PG (Procter & Gamble)
C (Citigroup)	INTC (Intel)	T (AT&T)
CSCO (Cisco)	JNJ (Johnson & Johnson)	UNH (United Health)
CVS (CVS Health)	JPM (JP Morgan)	UPS
CVX (Chevron)	K (Kellogg)	VZ (Verizon)
D (Dominion Energy)	MCK (McKesson)	WFC (Wells Fargo)
DIS (Disney)	MRK (Merck)	WMT (Walmart)
F (Ford)	MSFT (Microsoft)	XOM (Exxon)

There are 7 attributes: closing price (Close), market capitalisation in M$ (Mkt_Cap), price-to-book ratio (P2B), 1 month volatility (Vol_1M), debt-to-equity ratio (D2E) and profitability margin (Prof_Marg).
Finally, the time range is 2000-2021.

load("data_new.RData")                    # Loading the data: IF DIRECTORY OK!
data <- data %>% arrange(date,symbol)     # Ranked first according to date and then stocks
summary(data)                             # Descriptive statistics

##       date               symbol              close              vol_1y       
##  Min.   :2000-11-30   Length:8760        Min.   :  0.2525   Min.   :0.05741  
##  1st Qu.:2006-12-23   Class :character   1st Qu.: 29.3026   1st Qu.:0.15982  
##  Median :2013-01-15   Mode  :character   Median : 49.4400   Median :0.21048  
##  Mean   :2013-01-14                      Mean   : 75.0564   Mean   :0.24509  
##  3rd Qu.:2019-02-07                      3rd Qu.: 87.4050   3rd Qu.:0.29136  
##  Max.   :2025-02-28                      Max.   :960.0200   Max.   :1.59204  
##                                                                              
##     mkt_cap               p2b                d2e             revenue         
##  Min.   :5.024e+09   Min.   :  0.2767   Min.   :   0.00   Min.   :4.827e+09  
##  1st Qu.:6.452e+10   1st Qu.:  1.8193   1st Qu.:  35.59   1st Qu.:3.749e+10  
##  Median :1.374e+11   Median :  2.9471   Median :  70.67   Median :6.671e+10  
##  Mean   :1.855e+11   Mean   : 11.5858   Mean   : 248.47   Mean   :1.018e+11  
##  3rd Qu.:2.123e+11   3rd Qu.:  5.5292   3rd Qu.: 199.18   3rd Qu.:1.284e+11  
##  Max.   :3.463e+12   Max.   :540.0111   Max.   :5814.97   Max.   :6.810e+11  
##                      NA's   :121        NA's   :56        NA's   :1168       
##     scope_1         
##  Min.   :     9489  
##  1st Qu.:    95316  
##  Median :   567297  
##  Mean   :  9503221  
##  3rd Qu.:  1500000  
##  Max.   :145500000  
##  NA's   :2003

symbols <- unique(data$symbol)            # Set of assets
t_all <- unique(data$date)                # Set of dates
N <- n_distinct(data$symbol)              # Number of stocks

This simple table shows a possible outlier for the P2B variable. The maximum value is clearly out of range.

Next, we format the data for future use. Notably, we compute & store returns.

data <- data %>% 
  group_by(symbol) %>%                                      # Grouping: returns computed stock-by-stock
  mutate(return = close / dplyr::lag(close) - 1) %>%        # Adding returns
  ungroup()

returns <- data %>%                                           # Take data
  select(symbol, date, return) %>%                          # Select 3 columns
  pivot_wider(names_from = symbol, values_from = return) |> # Put them into 'matrix' format
  na.omit()
returns                                                       # Show the returns

## # A tibble: 291 × 31
##    date          AAPL      BA     BAC       C    CSCO      CVS      CVX        D
##    <date>       <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl>    <dbl>    <dbl>
##  1 2000-12-31 -0.0985 -0.0443  0.149   0.0251 -0.201   0.0538   0.0313   0.117  
##  2 2001-01-31  0.454  -0.114   0.173   0.0961 -0.0212 -0.0123  -0.0137  -0.0776 
##  3 2001-02-28 -0.156   0.0632 -0.0719 -0.121  -0.367   0.0304   0.0286   0.0608 
##  4 2001-03-31  0.209  -0.104   0.0961 -0.0854 -0.332  -0.0411   0.0250  -0.0166 
##  5 2001-04-30  0.155   0.109   0.0228  0.0927  0.0738  0.00786  0.0998   0.0624 
##  6 2001-05-31 -0.217   0.0176  0.0580  0.0427  0.134  -0.0687  -0.00528 -0.0320 
##  7 2001-06-30  0.165  -0.116   0.0132  0.0310 -0.0550 -0.297   -0.0578  -0.0931 
##  8 2001-07-31 -0.192   0.0527  0.0598 -0.0498  0.0560 -0.0671   0.00983  0.00599
##  9 2001-08-31 -0.0128 -0.125  -0.0333 -0.0888 -0.150   0.00278 -0.00700  0.0407 
## 10 2001-09-30 -0.164  -0.346  -0.0504 -0.115  -0.254  -0.0806  -0.0661  -0.0572 
## # ℹ 281 more rows
## # ℹ 22 more variables: DIS <dbl>, F <dbl>, FDX <dbl>, GE <dbl>, HD <dbl>,
## #   IBM <dbl>, INTC <dbl>, JNJ <dbl>, JPM <dbl>, K <dbl>, LLY <dbl>, MCK <dbl>,
## #   MRK <dbl>, MSFT <dbl>, ORCL <dbl>, PFE <dbl>, PG <dbl>, T <dbl>, VZ <dbl>,
## #   WFC <dbl>, WMT <dbl>, XOM <dbl>

3 What students often do…

3.1 Binary signal-based portfolios

Often, students think in terms of simple signals, e.g., based on crossings of moving averages (to detect trend changes):
- simple asset: this is complex because binary signals impose to buy/sell possibly frequently, which is costly
- multiple assets: even with several assets, binary signals are costly because positions are not very nuanced

Solution: real-valued signals. But adding a layer of optimization is preferable (include constraints!).

3.2 Static optimization

Now, let’s consider the usual Markowitz allocation subject to the budget constraint \(\textbf{w}'\textbf{1}=1\):

\[\text{max} \ \textbf{w}'\boldsymbol{\mu} - \frac{\gamma}{2}\textbf{w}'\boldsymbol{\Sigma}\textbf{w}, \quad \text{s.t.} \quad \textbf{w}'\textbf{1}=1 \] with Lagrangian: \[L(\textbf{w})=\textbf{w}'\boldsymbol{\mu} - \frac{\gamma}{2}\textbf{w}'\boldsymbol{\Sigma}\textbf{w}+\delta (\textbf{w}'\textbf{1}-1)\] so that \[\frac{\partial L}{\partial \textbf{w}}=\boldsymbol{\mu}-\gamma\boldsymbol{\Sigma}\textbf{w}+\delta\textbf{1}=0 \quad \Rightarrow \quad \textbf{w}=\gamma^{-1}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}+\delta \textbf{1}),\] where \(\delta\) is chosen to satisfy the budget constraint. A particular case is the minimum variance portfolio which corresponds to \(\boldsymbol{\mu}=\textbf{1}\): the optimal portfolio is proportional to \(\boldsymbol{\Sigma}^{-1} \textbf{1}\).

Sigma <- returns %>%                                       # Covariance matrix
  select(-date) %>% 
  cov(use = "complete.obs")
Sigma[1:9,1:9] %>% round(3)                                # A quick look at Sigma

##       AAPL    BA   BAC     C  CSCO   CVS   CVX     D   DIS
## AAPL 0.010 0.002 0.003 0.003 0.003 0.001 0.002 0.001 0.003
## BA   0.002 0.009 0.004 0.005 0.003 0.002 0.002 0.001 0.003
## BAC  0.003 0.004 0.012 0.010 0.003 0.002 0.003 0.001 0.004
## C    0.003 0.005 0.010 0.014 0.004 0.003 0.003 0.001 0.005
## CSCO 0.003 0.003 0.003 0.004 0.008 0.002 0.002 0.001 0.003
## CVS  0.001 0.002 0.002 0.003 0.002 0.006 0.002 0.001 0.002
## CVX  0.002 0.002 0.003 0.003 0.002 0.002 0.004 0.001 0.002
## D    0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.003 0.001
## DIS  0.003 0.003 0.004 0.005 0.003 0.002 0.002 0.001 0.006

Min_var_weights <- solve(Sigma) %*% rep(1,N)               # Proportional weights
Min_var_weights <- Min_var_weights / sum(Min_var_weights)  # Exact (scaled) weights
data.frame(symbols, Min_var_weights) %>%
  ggplot(aes(x = Min_var_weights, y = symbols)) +             # Plot
  geom_col(alpha = 0.5) + theme_light()

Naturally, the stocks with super low volatility (Disney, Walmart) are favored! Let’s see when the weights are driven by the mean returns… Mind the vector multiplication operator %*%.

mu <- returns %>%                                      # Mean vector 
  select(-date) %>%                                  # Remove date
  colMeans()
Marko_weights <- solve(Sigma) %*% mu                   # Proportional weights
Marko_weights <- Marko_weights / sum(Marko_weights)    # Exact (scaled) weights
bind_cols(symbols, Marko_weights) %>%
  ggplot(aes(x = Marko_weights, y = symbols)) +      # Plot
  geom_col(alpha = 0.5) + theme_light()

(Marko_weights[,1] %*% mu)/sqrt(t(Marko_weights) %*% Sigma %*% Marko_weights) # In-sample SR! 

##           [,1]
## [1,] 0.4854783

AAPL and UNH have high returns => high weights.

Now let’s look at 2 ways of processing these weights into portfolio values.

w_matrix <- matrix(Marko_weights, 
                   ncol = N, nrow = length(t_all) - 1, byrow = T) # Weight matrix
(rowSums((returns %>% select(-date)) * w_matrix) + 1) %>%         # Simple product = element-wise
  cumprod() %>% round(3)                                          # Cumulative product + rounding

##   [1]    1.090    1.209    1.162    1.266    1.366    1.376    1.393    1.424
##   [9]    1.403    1.346    1.317    1.418    1.475    1.568    1.579    1.719
##  [17]    1.776    1.695    1.660    1.522    1.638    1.617    1.531    1.446
##  [25]    1.492    1.497    1.428    1.445    1.460    1.603    1.714    1.682
##  [33]    1.765    1.671    1.754    1.632    1.718    1.663    1.799    1.921
##  [41]    1.978    2.083    2.256    2.257    2.418    2.641    2.941    3.328
##  [49]    3.388    3.749    3.929    3.718    3.545    3.644    3.727    4.095
##  [57]    4.391    4.628    4.976    5.417    5.689    5.932    5.747    5.780
##  [65]    5.904    5.656    5.635    5.802    5.838    6.389    6.607    6.799
##  [73]    6.591    6.897    6.888    7.209    7.469    7.923    7.616    7.902
##  [81]    8.116    8.639    9.751   10.377   10.724   10.789   10.239   11.756
##  [89]   11.964   12.857   12.409   13.118   13.292   12.704   12.581   11.763
##  [97]   12.559   13.319   13.432   12.335   12.927   13.599   14.673   15.477
## [105]   13.788   14.819   15.888   16.772   17.668   16.030   17.009   17.661
## [113]   18.340   17.837   17.836   17.479   16.001   17.576   17.919   17.255
## [121]   19.048   19.193   19.623   19.665   19.349   20.094   19.337   21.353
## [129]   22.194   20.777   20.756   20.475   21.507   22.355   24.899   26.283
## [137]   25.925   26.292   27.532   28.555   29.944   30.118   28.744   29.461
## [145]   27.494   27.431   27.813   27.677   27.263   28.260   26.167   28.907
## [153]   29.827   29.413   31.638   34.496   33.888   34.122   35.908   37.869
## [161]   38.318   40.405   40.445   39.530   41.891   42.918   45.489   46.956
## [169]   46.730   48.071   49.014   48.445   48.791   49.391   50.734   49.213
## [177]   44.995   44.733   44.249   45.986   45.790   44.970   42.694   45.478
## [185]   43.220   44.693   43.697   47.007   46.192   46.979   46.659   47.151
## [193]   49.308   53.619   58.771   60.225   60.380   63.139   64.990   68.021
## [201]   69.225   67.490   70.924   75.967   79.514   84.541   81.540   75.568
## [209]   75.378   78.561   79.689   86.996   96.871   98.335   97.062   96.827
## [217]   89.877   91.080   93.528   97.355  101.547   97.262  101.928  110.324
## [225]  119.288  120.611  119.027  125.815  125.391  134.445  128.801  126.954
## [233]  142.168  146.794  167.448  179.304  201.128  198.598  182.753  196.514
## [241]  199.309  226.986  220.354  216.024  239.029  237.609  259.219  284.890
## [249]  290.659  273.662  317.236  307.479  328.503  332.649  344.027  370.068
## [257]  383.618  352.209  353.748  395.330  407.911  405.275  463.101  484.198
## [265]  459.359  452.894  459.279  458.383  526.460  556.379  618.626  654.095
## [273]  694.044  664.224  693.502  722.200  693.480  736.293  810.436  776.098
## [281]  820.758  878.139  988.512  973.878 1132.218 1054.755 1088.193 1233.225
## [289] 1257.834 1215.720 1192.896

The above code is ugly. Because the data is neatly ordered, we can do better! Let’s add a new column to the data.

data %>% 
  select(symbol, date, return) %>%                          # Keep only necessary cols
  bind_cols(w_Marko = rep(Marko_weights, length(t_all)))    # Bind cols: beware to the order!

## # A tibble: 8,760 × 4
##    symbol date       return  w_Marko
##    <chr>  <date>      <dbl>    <dbl>
##  1 AAPL   2000-11-30     NA  0.374  
##  2 BA     2000-11-30     NA  0.0876 
##  3 BAC    2000-11-30     NA  0.0886 
##  4 C      2000-11-30     NA -0.306  
##  5 CSCO   2000-11-30     NA -0.0382 
##  6 CVS    2000-11-30     NA  0.0195 
##  7 CVX    2000-11-30     NA  0.136  
##  8 D      2000-11-30     NA -0.0404 
##  9 DIS    2000-11-30     NA  0.00359
## 10 F      2000-11-30     NA -0.0925 
## # ℹ 8,750 more rows

Now we can add benchmark weights and proceed.

p <- data %>% 
  select(symbol, date, return) %>%
  bind_cols(w_Marko = rep(Marko_weights, length(t_all))) %>%
  mutate(w_benchmark = 1/30) %>%                             # Add benchmark weights 1/N
  pivot_longer(cols = c(w_Marko, w_benchmark),               # Pivot to tidy
               names_to = "portfolio",  
               values_to = "weight") %>%
  group_by(date, portfolio) %>%                              # Group to compute portf return
  summarise(return = sum(weight * return)) %>%               # Portf return, date-by-date
  ungroup() %>%                                              # Ungroup
  na.omit() %>%                                              # Remove outliers
  group_by(portfolio) %>%                                    # Group along the 2 portfolios
  mutate(value = cumprod(1+return)) %>%                      # Compute portfolio values
  ggplot(aes(x = date, y = value, color = portfolio)) +      # Plot
  geom_line() + theme_light() + scale_y_log10() +
  scale_color_manual(values = c("#706C66", "#0A7BB0"), labels = c("Benchmark", "Markowitz")) + 
  theme(legend.position = c(0.3,0.7))

ggplotly(p, width = 900)

Given the log-scale: is this reasonably realistic? Obvious answer: NO! The weights have been chosen with knowledge of the trajectory!!! Also: portfolio managers update the weights…

What if we want to add constraints? Indeed, it’s often not possible to short weights (ex: with Citi…).
We will use a super simple, yet super useful package: quadprog!
quadprog solves the following programs:

\[\underset{\textbf{x}}{\min} \ -\textbf{d}'\textbf{x} + \frac{1}{2}\textbf{x}'\textbf{D}\textbf{x}=\underset{\textbf{x}}{\max} \ \textbf{d}'\textbf{x} - \frac{1}{2}\textbf{x}'\textbf{D}\textbf{x}, \quad \text{subject to} \ \textbf{A}'\textbf{x} \ge \textbf{b}_0,\] where in the constraints the inequality is meant element-wise and the first \(m\) inequalities are in fact equalities.
In our cases, it is clear that \(\textbf{d}\) is the vector of expected returns and \(\textbf{D}\) is the covariance matrix.
Moreover, we want to impose the budget constraint plus a non-negativity constraint for all assets, thus:

\[\textbf{A}'=\begin{bmatrix} 1 & 1 & \dots & 1 \\ 1 & 0 &\dots & 0 \\ 0 & 1 & \dots & 0 \\ 0 & 0 & \ddots &0 \\ 0 & 0 & \dots & 1 \end{bmatrix}, \quad \textbf{b}_0=\begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}\]

are the natural choices

if(!require(quadprog)){install.packages("quadprog")}
library(quadprog)
A <- rbind(rep(1, N),
           diag(N))
b_zero <- c(1, rep(0, N))
Marko_LO <- solve.QP(                               # Long-only Markowitz
  Dmat = Sigma,
  dvec = mu,
  Amat = t(A),
  bvec = b_zero,
  meq = 1         # only 1 equality = budget constraints.
)$solution
bind_cols(symbols, Marko_LO) %>%
  ggplot(aes(x = Marko_LO, y = symbols)) +             # Plot
  geom_col(alpha = 0.5) + theme_light()

(Marko_LO %*% mu)/sqrt(Marko_LO %*% Sigma %*% Marko_LO) # IS SR

##           [,1]
## [1,] 0.2814352

So much for diversification!
What if we impose diversification. For instance by setting a strictly positive weight for all stocks… In the example below, we ask for 1% of the portfolio at least in each firm.

b_zero <- c(1, rep(0.01, N))
Marko_LO <- solve.QP( # Long-only Markowitz
  Dmat = Sigma,
  dvec = mu,
  Amat = t(A),
  bvec = b_zero,
  meq = 1         # only 1 equality = budget constraints.
)$solution
bind_cols(symbols, Marko_LO) %>%
  ggplot(aes(x = Marko_LO, y = symbols)) +         # Plot
  geom_col(alpha = 0.5) + theme_light()

(Marko_LO %*% mu)/sqrt(Marko_LO %*% Sigma %*% Marko_LO) # IS SR

##          [,1]
## [1,] 0.282648

Apple crushes the competition.

4 A universe of functions

Backtesting portfolio strategies is usually simple at first, but becomes more intricate when many options are considered. To keep things simple, it is more efficient to work with functions. Functions help compartmentalise the different tasks of the process. We will build functions that compute portfolio weights and others that evaluate the performance metrics of the strategies. While we will usually work with a loop on backtesting dates, it is possible to build functions that directly generate portfolio returns (see below: the map() function).

By definition, a portfolio is a choice of weights that sum to one. Below, we implement two classical weighting schemes: the uniform portfolio (EW = equal weights) and the maximum Sharpe ratio portfolio (MSR). The latter is more complicated and requires inputs: namely the column vector of (expected) mean \(\mu\) and covariance matrix \(\Sigma\) of the assets. For simplicity, we will estimate them using sample moments - even though this is known to be a bad choice. This requires an additional argument in the function: the assets’ past returns. The MSR weights are \(w=\frac{\Sigma^{-1}\mu}{1'\Sigma^{-1}\mu}\). Both \(\mu\) and 1 are vectors here.

weights_msr <- function(returns){  # returns will refer to PAST returns
  m <- colMeans(returns)         # Vector of average returns
  sigma <- cov(returns)          # Covariance matrix
  w <- solve(sigma) %*% m        # Raw weights
  return(w / sum(w))             # Returns normalised weights
}

weights_ew <- function(returns){   # We keep the same syntax for simplicity
  N <- ncol(returns)             # Number of assets
  return(rep(1/N,N))             # Weight = 1/N
}

We are now ready to proceed with the initialisation of the variables that will use in the main loop.

sep_date <- as.Date("2010-01-01")                # This date separates in-sample vs out-of-sample
t_oos <- t_all[t_all > sep_date]                 # Out-of-sample dates (i.e., testing set)
portf_weights <- matrix(0, nrow = length(t_oos), # Will store portfolio weights
                        ncol = length(symbols)) 
portf_returns <- c()                             # Initialize portfolio returns

At last, we can proceed to the main loop. A note of caution: both in the weighting scheme functions and in the loop below, we assume well defined (finite, i.e., non NA) data. Obviously, feeding NA data in the system will produce NA outputs. There are only four steps in the loop:

1. extract the data
   
2. compute portfolio weights
   
3. compute realised returns
   
4. derive the return of the portfolio

for(t in 1:length(t_oos)){
  temp_data <- returns %>% filter(date < t_oos[t]) # 1) Extracting past data: expanding window! 
  portf_weights[t,] <- temp_data %>%               # 2) Take this past data 
    select(-date) %>%                                # Take out the date column
    na.omit() %>%                                    # Remove missing points
    weights_msr()                                    # Appply the weighting scheme function
  realised_returns <- returns %>%                  # 3) Take returns
    filter(date ==  t_oos[t]) %>%                    # Keep only current date
    select(-date)                                    # Take out the date column
  portf_returns[t] <- sum(portf_weights[t,] * realised_returns) # 4) Compute the return
}

In the above backtest, the amount of data considered to form the portfolio decision increases with time (expanding window). It is easy to fix the number of data point at each step and proceed on rolling windows (exercise). Likewise, switching from MSR to EW is immediate (just change the weighting function).

The main output is the vector of portfolio returns. We can easily plot the evolution of the portfolio through time.

port <- data.frame(date = t_oos, 
                   portfolio = cumprod(1+portf_returns))    # Portf. values via cumulative product
port %>% ggplot() + theme_light() +
  geom_line(aes(x = date, y = portfolio))                   # Plot

Below, we create a function that takes a series of returns as input and provides a few simple performance metrics.

perf_met <- function(returns){
  avg_ret <- mean(returns, na.rm = T)*12                  # Arithmetic mean 
  vol <- sd(returns, na.rm = T)*sqrt(12)                  # Volatility
  Sharpe_ratio <- avg_ret / vol                           # Sharpe ratio
  VaR_5 <- quantile(returns, 0.05)                        # Value-at-risk
  met <- data.frame(avg_ret, vol, Sharpe_ratio, VaR_5)    # Aggregation of all of this
  rownames(met) <- "metrics"
  return(met)
}
perf_met(portf_returns) # Let's test on the actual returns

##           avg_ret       vol Sharpe_ratio       VaR_5
## metrics 0.1774365 0.2175875    0.8154721 -0.07412241

Note: the values are not annualised. The annualisation can be directly coded into the function. A heuristic way to proceed is to multiply average returns by 12 and volatilities by \(\sqrt{12}\) - though this omits the compounding effect. These simplifications are ok if they are used to compare strategies. The Value at Risk is nonetheless dependent on the return horizon and cannot be proxied so simply.

An important indicator is the turnover of the portfolio: it assesses asset rotation and thus impacts transaction costs. Simple turnover computes average absolute variation in portfolio weights: \(\text{Turn}=\frac{1}{T}\sum_{t=2}^T\sum_{n=1}^N|w_t^n-w_{t-1}^n|\). Full turnover takes into account the variation of the weights between rebalancing dates: \(\text{Turn}=\frac{1}{T}\sum_{t=2}^T\sum_{n=1}^N|w_t^n-w_{t-}^n|\), where \(t-\) is the time just before the rebalancing date.

turnover_simple <- function(weights, t_oos){
  turn <- 0
  for(t in 2:length(t_oos)){                             # Loop on dates
    turn <- turn + sum(abs(weights[t,] - weights[t-1,])) # Variation in weights
  }
  return(12*turn/(length(t_oos)-1))                         # BEWARE: monthly value!
}
turnover_simple(portf_weights, t_oos)                    # Simple turnover!

## [1] 3.576248

Below, we switch to the full definition of the turnover. In this case, the variation in weights is adjusted because the weight prior to rebalancing have drifted.

turnover_full <- function(weights, asset_returns, t_oos){
  turn <- 0
  for(t in 2:length(t_oos)){
    realised_returns <- returns %>% filter(date == t_oos[t]) %>% select(-date)
    prior_weights <- weights[t-1,] * (1 + realised_returns)       # Before rebalancing
    turn <- turn + apply(abs(weights[t,] - prior_weights/sum(prior_weights)),1,sum)
  }
  return(12*turn/(length(t_oos)-1))
}
asset_returns <- returns %>% filter(date > sep_date)    # Asset returns over the experiment
turnover_full(portf_weights, asset_returns, t_oos)      # Real turnover

## [1] 5.282898

Note that the turnover is computed at the monthly frequency. The second value is always more realistic; in this case it is substantially higher compared to the simplified proxy. For the sake of compactness, the turnover should be included in the perf_met() function.

5 Extensions

5.1 Comparing strategies

An important generalisation of the above framework is to consider more than one strategy. Comparisons are commonplace in the asset management industry: obviously, people look for the best strategy (according to particular goals, beliefs, and preferences). Below, we show how this can be handled. We will compare the two strategies that we mentionned above. First, we need to re-initiate the variables because their dimension will change (NOTE: an alternative route would be to work with lists).

Tt <- length(t_oos)                                              # Nb of computation dates 
# Avoid T because T = TRUE!
nb_port <- 2                                                     # Nb of portfolios
portf_weights <- array(0, dim = c(Tt, nb_port, length(symbols))) # Store weights
portf_returns <- matrix(0, nrow = Tt, ncol = nb_port)            # Store returns

Second, we embed all weighting schemes into one single function.

weights_multi <- function(returns, j){   # Strategies are indexed by j
  if(j == 1){ # j = 1 => MSR
    return(weights_msr(returns))
  }
  if(j == 2){ # j = 2 => EW
    N <- length(returns[1,])
    return(rep(1/N,N))
  }
}

We decided to recode the EW strategy, but we could have used the weights_ew() function instead. Finally, the main loop is only marginally different from the single strategy loop.

for(t in 1:length(t_oos)){
  temp_data <- returns %>% 
    filter(date < t_oos[t]) %>%                         # Same first step
    select(-date)
  realised_returns <- returns %>%                       # Same third step: returns
    filter(date ==  t_oos[t]) %>% 
    select(-date)
  for(j in 1:nb_port){                                  # This is the novelty: we loop on the strategies 
    portf_weights[t,j,] <- weights_multi(temp_data, j)  # The weights' function is indexed by j
    portf_returns[t,j] <- sum(portf_weights[t,j,] * realised_returns)
  }
}

apply(portf_returns, 2, perf_met) %>%                   # Taking perf metrics
  unlist() %>%                                          # Flattening the list
  matrix(nrow = 2, byrow = T) %>%                       # Ordering them
  `colnames<-`(c("avg_ret", "vol", "SR", "VaR")) %>%    # Adding column names
  data.frame()                                          # Converting to dataframe

##     avg_ret       vol        SR         VaR
## 1 0.1774365 0.2175875 0.8154721 -0.07412241
## 2 0.1091758 0.1358828 0.8034553 -0.05443078

apply(portf_weights, 2, turnover_simple, t_oos = t_oos) %>% unlist()

## [1] 3.576248 0.000000

We recall the order of strategies: MSR first (line) and EW second (line).

Again, turnover should (and will) be added to the performance metric function. Since uniform weights are constant, their simplified turnover is zero. In practice, that is not the case because weights evolve according to asset returns. The adjustment would imply only a small turnover.

There is a well-documented substantial difference between the two weigting schemes in terms of asset rotation. The MSR that we compute is clearly not competitive: even before transaction costs, its Sharpe ratio is smaller than that of the EW portfolio (see DeMiguel et al. (2009) for further evidence on the robustness of the EW portfolio).

5.2 Make do without loops: functional programming

Finally, we show how to bypass the loop over dates. While this is not useful on small datasets, it can save time on large databases because loops are notoriously time-consuming. Below, we show how to proceed with the map() function in the simple case with only one strategy.

# We create a function that will compute returns for each date:
port_map <- function(t_oos, returns, weight_func){ # Weighting function as argument!             
  temp_data <- returns %>% filter(date < t_oos)    # Still expanding window...
  portf_weights <- temp_data %>% 
    select(-date) %>% 
    weight_func()
  realised_returns <- returns %>% 
    filter(date ==  t_oos) %>% 
    select(-date)
  return(sum(portf_weights * realised_returns))
}
# the map() function does it all!
portf_returns <- t_oos %>%                         # The variable over which we loop
  map(~port_map(.x, 
                returns = returns, 
                weight_func = weights_msr)         # MSR weights
  ) %>%                                            # Is sent to the map() function
  unlist()                                         # The output is flattened
perf_met(portf_returns)                            # Compute the perf metrics

##           avg_ret       vol Sharpe_ratio       VaR_5
## metrics 0.1774365 0.2175875    0.8154721 -0.07412241

6 Characteristics-based portfolio choice

In this section, we start the chapter of strategies based on firm characteristics (features). As an illustration, we check a well-documented (though still controversial) anomaly: the size effect. We build two portfolios: in the first (resp. second) one, we invest in the firms that have a below (resp. above) median market capitalisation. The first portfolio will be a ‘small’ portfolio and the second one a ‘large’ one. Stocks are equally-weighted inside the portfolios.

First, we prepare the variables and define the weight function.

nb_port <- 2                                                         # Nb of portfolios
portf_weights <- array(0, dim = c(Tt-1, nb_port, length(symbols)))   # Where we store weights
portf_returns <- matrix(0, nrow = Tt-1, ncol = nb_port)              # Where we store returns

weights_cap <- function(data,j){    
  # More general than before: we feed all the data, not just returns
  m <- median(data$mkt_cap)       # Compute the median market cap
  n <- nrow(data)                 # Compute the number of assets
  if(j == 1){return((data$mkt_cap < m)/n*2)} # Small cap
  if(j == 2){return((data$mkt_cap > m)/n*2)} # Large cap
}

There will only be Tt-1 dates since we lose one because of the computation of future returns.

Second, we can launch the backtesting loop.

for(t in 2:length(t_oos)){
  temp_data <- data %>% filter(date == t_oos[t-1])        # We keep the data of the previous date
  for(j in 1:nb_port){                                    # We loop on the strategies 
    portf_weights[t-1,j,] <- weights_cap(temp_data, j)    # The weights' function is indexed by j
    realised_returns <- returns %>% 
      filter(date ==  t_oos[t]) %>%                       # Shift in the date
      select(-date)
    portf_returns[t-1,j] <- sum(portf_weights[t-1,j,] * realised_returns)
  }
}

Third, we proceed to performance metrics.

apply(portf_returns, 2, perf_met) %>%                   # Taking perf metrics
  unlist() %>%                                          # Flattening the list
  matrix(nrow = 2, byrow = T) %>%                       # Ordering them
  `colnames<-`(c("avg_ret", "vol", "SR", "VaR")) %>%    # Adding column names
  data.frame()                                          # Converting to dataframe

##      avg_ret       vol        SR         VaR
## 1 0.13142480 0.1464983 0.8971082 -0.05777850
## 2 0.09045546 0.1373620 0.6585188 -0.05916399

Small firms do indeed generate a higher level of performance! In order to reach this conclusion in a rigourous fashion, we would need to perform the same analysis on at least 1,000 stocks (ideally, more) and on 5 to 10 portfolio sorts (from very small firms to very large ones). Below, we check the weights of the portfolio on one particular date.

small <- portf_weights[2,1,] # t = 2, j = 1 (t_oos[2] = 2010-03-01, small firms)
large <- portf_weights[2,2,] # t = 2, j = 2 (t_oos[2] = 2010-03-01, large firms)
data.frame(small, large, row.names = symbols)

##           small      large
## AAPL 0.00000000 0.06666667
## BA   0.06666667 0.00000000
## BAC  0.00000000 0.06666667
## C    0.06666667 0.00000000
## CSCO 0.00000000 0.06666667
## CVS  0.06666667 0.00000000
## CVX  0.00000000 0.06666667
## D    0.06666667 0.00000000
## DIS  0.06666667 0.00000000
## F    0.06666667 0.00000000
## FDX  0.06666667 0.00000000
## GE   0.06666667 0.00000000
## HD   0.06666667 0.00000000
## IBM  0.00000000 0.06666667
## INTC 0.00000000 0.06666667
## JNJ  0.00000000 0.06666667
## JPM  0.00000000 0.06666667
## K    0.06666667 0.00000000
## LLY  0.06666667 0.00000000
## MCK  0.06666667 0.00000000
## MRK  0.06666667 0.00000000
## MSFT 0.00000000 0.06666667
## ORCL 0.06666667 0.00000000
## PFE  0.00000000 0.06666667
## PG   0.00000000 0.06666667
## T    0.00000000 0.06666667
## VZ   0.06666667 0.00000000
## WFC  0.00000000 0.06666667
## WMT  0.00000000 0.06666667
## XOM  0.00000000 0.06666667

Indeed, some stocks have zero weights and others 1/15.

Finally, let’s see how we could have coded those strategies using the tidyverse (and a lot of piping!).

data %>% filter(date > sep_date) %>%                        # Keep only the out-of-sample backtesting dates
  group_by(symbol) %>%                                    # Group by stock
  mutate(f_return = dplyr::lead(return)) %>%              # Compute forward (i.e., realised) return
  na.omit() %>%                                           # Take out NAs
  group_by(date) %>%                                      # Group by dates
  mutate(mkt_cap_binary = mkt_cap < median(mkt_cap)) %>%  # Compute median cap for each date
  ungroup() %>%
  group_by(mkt_cap_binary) %>%                            # Group by Mkt_Cap: small vs large
  summarise(avg_return = mean(f_return),
            vol = sd(f_return))                 # Simple pivot table

## # A tibble: 2 × 3
##   mkt_cap_binary avg_return    vol
##   <lgl>               <dbl>  <dbl>
## 1 FALSE             0.00795 0.0628
## 2 TRUE              0.0102  0.0738

It can be useful to see how often stocks switch from one family to another (from below median to above median or vice-versa). Below, we show a plot of Mkt_Cap, conditional on Mkt_Cap being above the current median.

data %>% group_by(date) %>%                              # Group by date
  mutate(mkt_cap_median = median(mkt_cap)) %>%         # Compute median cap for each date
  filter(mkt_cap > mkt_cap_median) %>%                 # Keep only the large stocks
  ggplot(aes(x = date, y = mkt_cap, color = symbol)) + # Plot
  geom_line() + scale_y_log10() +
  geom_line(aes(x = date, y = mkt_cap_median), color = "black") 

# The black line shows the running median

The straight lines show the discontinuities: one stock being large at some point in time, then small and then large again. The straight lines show the periods when the stock was small. The black line shows the median capitalisation (in the sample). Finally, because we focus in the zone close to the median and impose an upper limit of 250B$, there are some missing points.

7 Sustainable investing

ESG investing has gained a lot of traction (see my review paper: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3715753).
In this section, we briefly explore this topic.

First, a glance at the data. Let’s look at the distribution of CO2 intensity (from Sustainalytics).
We must use a log-scale because the differences are pronounced.

data %>%
  mutate(year = year(date),
         co2_int = scope_1/revenue) %>%
  ggplot(aes(x = co2_int, fill = as.factor(year))) + 
  geom_histogram() + labs(fill = "year") + theme_light() +
  scale_x_log10()

Let’s continue with a “simple” piping flow that creates a compact dataset.

data_esg <- data %>%
  mutate(year = year(date),
         co2_int = scope_1/revenue) |>
  na.omit() |>
  filter(date > "2010-02-01") %>%                    # Data available after 2014-02-01
  group_by(date) |>
  mutate(int_score = ecdf(co2_int)(co2_int),
         ESG_group = if_else(int_score > median(int_score, na.rm=T), "brown", "green")) |>
  ungroup() |>
  group_by(symbol) %>%
  mutate(return = dplyr::lead(close)/close - 1) %>%  # Forward return
  ungroup() 
data_esg

## # A tibble: 4,696 × 14
##    date       symbol close vol_1y   mkt_cap    p2b    d2e revenue scope_1 return
##    <date>     <chr>  <dbl>  <dbl>     <dbl>  <dbl>  <dbl>   <dbl>   <dbl>  <dbl>
##  1 2010-02-28 AAPL    7.31  0.266   1.63e11  5.19     0   4.29e10  2.45e4 0.148 
##  2 2010-02-28 BA     63.2   0.288   3.93e10 10.9    581.  6.83e10  5.79e5 0.150 
##  3 2010-02-28 CSCO   24.3   0.287   1.26e11  3.28    26.6 3.61e10  5.26e4 0.0699
##  4 2010-02-28 CVS    33.8   0.271   4.54e10  1.25    31.2 9.82e10  1.57e5 0.0833
##  5 2010-02-28 CVX    72.3   0.166   1.54e11  1.68    11.4 1.67e11  6.08e7 0.0488
##  6 2010-02-28 D      38.0   0.131   2.32e10  2.08   158.  1.44e10  5.90e7 0.0821
##  7 2010-02-28 DIS    31.2   0.257   5.06e10  1.47    36.5 3.61e10  5.43e5 0.117 
##  8 2010-02-28 F      10.5   1.28    3.04e10  2.44  4788.  1.16e11  1.5 e6 0.0707
##  9 2010-02-28 FDX    84.8   0.356   1.73e10  1.27    20.1 3.55e10  1.41e7 0.102 
## 10 2010-02-28 GE     76.9   0.430   9.64e10  0.823  402.  1.54e11  2.70e6 0.133 
## # ℹ 4,686 more rows
## # ℹ 4 more variables: year <dbl>, co2_int <dbl>, int_score <dbl>,
## #   ESG_group <chr>

Based on this short dataset, let’s have a quick look at the perf of the 2 groups (high ESG versus low ESG).

data_esg %>%
  na.omit() |>
  group_by(date, ESG_group) %>%
  summarise(return = mean(return, na.rm = T)) %>%
  group_by(ESG_group) %>%
  summarise(avg_return = 12*mean(return, na.rm = T),  # Be careful with na.rm = T
            vol = sqrt(12)*sd(return, na.rm = T),
            SR = avg_return/vol)

## # A tibble: 2 × 4
##   ESG_group avg_return   vol    SR
##   <chr>          <dbl> <dbl> <dbl>
## 1 brown         0.0877 0.131 0.671
## 2 green         0.128  0.135 0.947

Two conclusions:
1. High ESG does slightly better.
2. But the difference is not huge (t-stat => exercise!)

From individual annual returns: is there are pattern?

data_esg %>%
  group_by(year, symbol) %>%
  summarise(avg_int = mean(int_score, na.rm = T),
            avg_return = mean(return, na.rm = T)) %>%
  ggplot(aes(x = avg_int, y = avg_return)) + geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") + theme_light()

No pattern at all, it seems…

Let’s dig deeper and have a look at screening/selection intensity. We compute the average returns of firms that have ESG_rank above a given threshold. To do that, let’s make a function.

library(scales) # This will be used in the plot below to zoom on the values with special limits
ESG_impact <- function(level, data_esg){
  data_esg %>%
    filter(int_score > level) %>%   # Filter to keep only the firms with ESG > level
    pull(return) %>%
    mean(na.rm = T)
}

level_int <- seq(0.1, 0.9, by = 0.1)                # ESG_rank threshold
map_dbl(level_int, ESG_impact, data_esg = data_esg) # The raw result

## [1] 0.008228179 0.007404221 0.007398298 0.007602926 0.007295227 0.007692715
## [7] 0.007426928 0.006044294 0.005379577

level_int %>%
  bind_cols(return = map_dbl(level_int, ESG_impact, data_esg = data_esg)) %>%
  ggplot(aes(x = level_int, y = return)) + geom_col(alpha = 0.5) + 
  theme_light() #+ scale_y_continuous(limits=c(0.008,0.012), oob = rescale_none)

Overall, lower intensity selection seems promising indeed (but let’s not jump to general conclusions). Finally, let’s have a look at more homogeneous groups.

ESG_quant <- function(level, width, data_esg){ # A new function
  data_esg %>%
    filter(int_score > level - width/2,
           int_score < level + width/2) %>%     # Filter to keep firms with ESG around level
    pull(return) %>%
    mean(na.rm = T)
}
level_esg <- quantile(data_esg$int_score, seq(0.1,0.9, by = 0.1))
level_esg %>%
  bind_cols(return = map_dbl(level_esg, ESG_quant, width = 0.1, data_esg = data_esg)) %>%
  ggplot(aes(x = level_esg, y = return)) + geom_col(alpha = 0.5) + 
  theme_light() #+ scale_y_continuous(limits=c(0.007,0.0125), oob = rescale_none)

## New names:
## • `` -> `...1`

The effect is not monotonic: things are not that simple…

8 Exercises

8.1 Rolling window

Change the main loop so that only 60 points of data are given to the weighting scheme(s). Sixty points amount to 5 years of monthly data.

8.2 Other performance metrics

Using the PerformanceAnalytics package (install it first!), compute the maximum drawdown via: https://www.rdocumentation.org/packages/PerformanceAnalytics/versions/1.5.2/topics/maxDrawdown

Compute the transaction cost-adjusted SR: \(TC-SR=(\bar{r}-0.005*Turn)/\sigma\).

Add both metrics to the perf_met() function.

8.3 Minimum variance

Add the MV portfolio to the set of strategies. The weights depend only on the covariance matrix: \(w=\frac{\Sigma^{-1}1}{1'\Sigma^{-1}1}\).

8.4 map() expertise

Extend the map() syntax to the case with many strategies. BONUS: have a look at pmap()!

8.5 Realistic portfolios

In practice, many saveguards are applied, if only to reduce turnover. One such example is box constraint: the weights in the portfolio must not lie above or below user-specified thresholds.

In minimum variance, it is possible to reduce leverage by shrinking the covariance matrix towards the identity functions: \[\Sigma^*=(1-\alpha) \Sigma^S+\alpha I,\] where \(\Sigma^S\) is the sample covariance matrix and \(I\) the identity matrix. As \(\alpha\) increases to 1, the weights converge to the equally-weighted portfolio.