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 <- "2022-05-05"
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 <- 63                                       # 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), 1 month relative strength index (RSI_1M), debt-to-equity ratio (D2E) and profitability margin (Prof_Marg).
Finally, the time range is 2000-2021.

load("data.RData")                      # Loading the data: IF DIRECTORY OK!
data <- data %>% arrange(Date,Tick)     # Ranked first according to date and then stocks
summary(data)                           # Descriptive statistics
##       Tick           Date                Close             Vol_1M       
##  AAPL   : 255   Min.   :1999-12-31   Min.   :  0.217   Min.   :  5.867  
##  BA     : 255   1st Qu.:2005-03-31   1st Qu.: 19.166   1st Qu.: 15.970  
##  BAC    : 255   Median :2010-07-30   Median : 32.307   Median : 21.835  
##  C      : 255   Mean   :2010-07-30   Mean   : 49.261   Mean   : 27.054  
##  CSCO   : 255   3rd Qu.:2015-11-30   3rd Qu.: 57.307   3rd Qu.: 31.894  
##  CVS    : 255   Max.   :2021-02-26   Max.   :442.987   Max.   :268.563  
##  (Other):6120                                                           
##     Mkt_Cap             P2B                 D2E             Prof_Marg       
##  Min.   :   4467   Min.   :   0.1013   Min.   :    0.00   Min.   :-106.260  
##  1st Qu.:  67392   1st Qu.:   1.4223   1st Qu.:   31.83   1st Qu.:   5.008  
##  Median : 137516   Median :   2.3494   Median :   66.76   Median :  10.312  
##  Mean   : 160615   Mean   :  11.1931   Mean   :  218.70   Mean   :  11.712  
##  3rd Qu.: 210966   3rd Qu.:   4.2259   3rd Qu.:  213.09   3rd Qu.:  19.171  
##  Max.   :2255969   Max.   :1678.9053   Max.   :14585.61   Max.   : 145.583  
##                                                                             
##     ESG_rank      
##  Min.   :  7.317  
##  1st Qu.: 54.237  
##  Median : 72.222  
##  Mean   : 68.352  
##  3rd Qu.: 86.364  
##  Max.   :100.000  
##  NA's   :5214
tick <- unique(data$Tick)               # Set of assets
t_all <- unique(data$Date)              # Set of dates
N <- n_distinct(data$Tick)              # 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(Tick) %>%                                   # Grouping: returns computed stock-by-stock
    mutate(Return = Close / dplyr::lag(Close) - 1) %>%   # Adding returns
    ungroup()

returns <- data %>%                                    # Take data
    select(Tick, Date, Return) %>%                       # Select 3 columns
    pivot_wider(names_from = Tick, values_from = Return) # Put them into 'matrix' format
returns                                                # Show the returns
## # A tibble: 255 Ɨ 31
##    Date           AAPL      BA     BAC       C    CSCO      CVS      CVX       D
##    <date>        <dbl>   <dbl>   <dbl>   <dbl>   <dbl>    <dbl>    <dbl>   <dbl>
##  1 1999-12-31 NA       NA      NA      NA      NA      NA       NA       NA     
##  2 2000-01-31  0.00885  0.0739 -0.0349  0.0236  0.0222 -0.122   -0.0346   0.0637
##  3 2000-02-29  0.105   -0.167  -0.0503 -0.0896  0.207   0.00179 -0.0993  -0.106 
##  4 2000-03-31  0.184    0.0237  0.152   0.157   0.170   0.0732   0.238    0.0477
##  5 2000-04-28 -0.0862   0.0496 -0.0656 -0.0121 -0.103   0.160   -0.0791   0.171 
##  6 2000-05-31 -0.323   -0.0121  0.142   0.0540 -0.179   0        0.0936   0.0313
##  7 2000-06-30  0.248    0.0704 -0.224  -0.0312  0.116  -0.0805  -0.0825  -0.0628
##  8 2000-07-31 -0.0298   0.167   0.102   0.170   0.0295 -0.0112  -0.0687   0.0598
##  9 2000-08-31  0.198    0.102   0.142   0.107   0.0487 -0.0596   0.0785   0.179 
## 10 2000-09-29 -0.577    0.202  -0.0225 -0.0742 -0.195   0.247    0.00864  0.0982
## # ā€¦ with 245 more rows, and 22 more variables: DIS <dbl>, F <dbl>, GE <dbl>,
## #   HD <dbl>, IBM <dbl>, INTC <dbl>, JNJ <dbl>, JPM <dbl>, K <dbl>, MCK <dbl>,
## #   MRK <dbl>, MSFT <dbl>, ORCL <dbl>, PFE <dbl>, PG <dbl>, T <dbl>, UNH <dbl>,
## #   UPS <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.013 0.002 0.003 0.004 0.004 0.001 0.002 0.000 0.003
## BA   0.002 0.009 0.004 0.006 0.003 0.002 0.003 0.001 0.003
## BAC  0.003 0.004 0.013 0.011 0.003 0.002 0.002 0.001 0.004
## C    0.004 0.006 0.011 0.015 0.005 0.002 0.003 0.001 0.005
## CSCO 0.004 0.003 0.003 0.005 0.009 0.001 0.002 0.000 0.003
## CVS  0.001 0.002 0.002 0.002 0.001 0.006 0.002 0.001 0.002
## CVX  0.002 0.003 0.002 0.003 0.002 0.002 0.004 0.001 0.002
## D    0.000 0.001 0.001 0.001 0.000 0.001 0.001 0.002 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
bind_cols(tick, Min_var_weights) %>%
    ggplot(aes(x = Min_var_weights, y = tick)) +             # 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ā€¦

mu <- returns %>%                                      # Mean vector 
    select(-Date) %>%                                  # Remove date
    colMeans(na.rm = T)
Marko_weights <- solve(Sigma) %*% mu                   # Proportional weights
Marko_weights <- Marko_weights / sum(Marko_weights)    # Exact (scaled) weights
bind_cols(tick, Marko_weights) %>%
    ggplot(aes(x = Marko_weights, y = tick)) +         # 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.4745629

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
    na.omit() %>%                                                   # Remove missing points
    cumprod() %>% round(3)                                          # Cumulative product + rounding
##   [1]   0.957   0.895   0.931   0.933   0.927   0.996   0.985   1.097   1.085
##  [10]   1.095   1.101   1.209   1.236   1.256   1.300   1.359   1.310   1.316
##  [19]   1.412   1.411   1.312   1.274   1.359   1.405   1.516   1.526   1.673
##  [28]   1.898   1.879   1.858   1.702   1.729   1.671   1.593   1.538   1.607
##  [37]   1.639   1.564   1.597   1.630   1.781   1.846   1.871   1.886   1.835
##  [46]   1.928   1.859   1.954   1.960   2.043   2.138   2.111   2.240   2.366
##  [55]   2.368   2.516   2.706   2.904   3.293   3.341   3.570   3.697   3.769
##  [64]   3.665   3.762   3.857   4.055   4.219   4.486   4.718   4.957   5.056
##  [73]   5.044   5.060   4.884   4.965   4.759   4.843   5.204   5.171   5.360
##  [82]   5.560   5.720   5.700   5.829   5.894   6.004   6.103   6.269   6.035
##  [91]   6.066   6.167   6.389   7.167   8.004   8.508   7.763   7.271   7.284
## [100]   7.540   8.156   7.465   7.650   8.099   7.535   7.358   7.195   7.863
## [109]   8.345   7.881   7.529   7.940   8.635   9.435   9.963   8.916   9.034
## [118]   9.398  10.201  11.351  11.002  11.670  12.086  11.940  11.344  11.201
## [127]  11.591  11.362  12.097  12.296  11.989  12.214  12.611  13.273  13.464
## [136]  13.812  14.252  14.058  14.536  14.963  14.356  14.518  14.638  15.176
## [145]  15.599  16.918  17.804  17.594  18.055  19.206  18.919  19.155  18.661
## [154]  17.938  18.058  17.509  17.361  17.620  18.328  19.586  19.818  20.029
## [163]  21.619  22.017  22.751  23.817  25.805  25.872  26.049  27.052  29.119
## [172]  28.340  29.926  30.933  29.912  31.684  32.336  34.675  36.337  37.271
## [181]  37.692  38.300  37.915  37.414  38.274  38.052  38.422  35.745  36.012
## [190]  37.135  36.997  37.023  37.860  37.233  39.859  38.358  39.147  40.790
## [199]  42.238  41.284  42.089  43.308  45.886  47.406  50.376  53.211  53.917
## [208]  56.048  58.458  59.210  63.626  65.449  63.280  69.919  76.063  78.350
## [217]  83.468  80.938  76.859  77.800  79.071  80.365  86.609  92.206  94.659
## [226]  96.311 102.447  90.073  87.358  92.111  95.476  99.186  95.228  99.665
## [235]  97.140 103.540 104.853 108.992 110.806 115.953 111.383  97.044  89.036
## [244] 102.187 103.793  99.388  95.844 100.039  93.242  81.467  99.422  99.446
## [253]  99.910  98.918

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(Tick, Date, Return) %>%                          # Keep only necessary cols
    bind_cols(w_Marko = rep(Marko_weights, length(t_all)))  # Bind cols: beware to the order!
## # A tibble: 7,650 Ɨ 4
##    Tick  Date       Return w_Marko
##    <fct> <date>      <dbl>   <dbl>
##  1 AAPL  1999-12-31     NA  0.187 
##  2 BA    1999-12-31     NA  0.114 
##  3 BAC   1999-12-31     NA  0.117 
##  4 C     1999-12-31     NA -0.312 
##  5 CSCO  1999-12-31     NA -0.0181
##  6 CVS   1999-12-31     NA -0.0139
##  7 CVX   1999-12-31     NA  0.0668
##  8 D     1999-12-31     NA  0.228 
##  9 DIS   1999-12-31     NA  0.124 
## 10 F     1999-12-31     NA -0.0509
## # ā€¦ with 7,640 more rows

Now we can add benchmark weights and proceed.

p <- data %>% 
    select(Tick, 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(tick, Marko_LO) %>%
  ggplot(aes(x = Marko_LO, y = tick)) +             # Plot
  geom_col(alpha = 0.5) + theme_light()

(Marko_LO %*% mu)/sqrt(Marko_LO %*% Sigma %*% Marko_LO) # IS SR
##           [,1]
## [1,] 0.2714241

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(tick, Marko_LO) %>%
  ggplot(aes(x = Marko_LO, y = tick)) +         # Plot
  geom_col(alpha = 0.5) + theme_light()

(Marko_LO %*% mu)/sqrt(Marko_LO %*% Sigma %*% Marko_LO) # IS SR
##           [,1]
## [1,] 0.2495396

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 <- apply(returns, 2, mean)     # 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 <- length(returns[1,])         # 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(tick)) 
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)                     # Arithmetic mean 
    vol <- sd(returns, na.rm = T)                           # 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.010521 0.04392896    0.2395004 -0.06412108

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(turn/(length(t_oos)-1))                         # BEWARE: monthly value!
}
turnover_simple(portf_weights, t_oos)                    # Simple turnover!
## [1] 0.2189094

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(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] 0.3014361

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(tick)))   # 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
        na.omit() %>% 
        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.01052100 0.04392896 0.2395004 -0.06412108
## 2 0.01099715 0.03975785 0.2766032 -0.05482998
apply(portf_weights, 2, turnover_simple, t_oos = t_oos) %>% unlist()
## [1] 0.2189094 0.0000000

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) %>% 
        na.omit() %>%
        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.010521 0.04392896    0.2395004 -0.06412108

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(tick)))   # 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
            na.omit() %>% 
            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.013091631 0.04412177 0.2967159 -0.05881253
## 2 0.009295396 0.03820161 0.2433247 -0.05804917

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 = tick)
##           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
## GE   0.00000000 0.06666667
## HD   0.06666667 0.00000000
## IBM  0.00000000 0.06666667
## INTC 0.06666667 0.00000000
## JNJ  0.00000000 0.06666667
## JPM  0.00000000 0.06666667
## K    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
## UNH  0.06666667 0.00000000
## UPS  0.06666667 0.00000000
## 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(Tick) %>%                                      # 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))                  # Simple pivot table
## # A tibble: 2 Ɨ 2
##   Mkt_Cap_Binary avg_return
##   <lgl>               <dbl>
## 1 FALSE             0.00812
## 2 TRUE              0.0104

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 = Tick)) + # Plot
    geom_line() + ylim(75000,250000) +
    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 ESG_rank (from Sustainalytics).

load("data_esg.RData") # Loading the data
data_esg %>%
    mutate(Year = year(Date)) %>%
    ggplot(aes(x = ESG_rank, fill = as.factor(Year))) + 
    geom_histogram() + labs(fill = "Year") + theme_light()

There are few firms with extreme ESG_rank (below 10 or above 90).

Letā€™s continue with a ā€œsimpleā€ piping flow that creates a compact dataset.

data_esg2 <- data_esg %>%
    select(Tick, Date, Close, ESG_rank) %>%            # Keep relevant columns
    filter(Date > "2014-02-01") %>%                    # Data available after 2014-02-01
    group_by(Tick) %>%
    mutate(Return = dplyr::lead(Close)/Close - 1) %>%  # Forward return
    ungroup() %>%
    group_by(Date) %>%                                 # For each date: create ESG groups
    mutate(ESG_group = if_else(ESG_rank > median(ESG_rank), "High_ESG", "Low_ESG")) %>%
    ungroup()
data_esg2
## # A tibble: 32,640 Ɨ 6
##    Tick  Date       Close ESG_rank    Return ESG_group
##    <chr> <date>     <dbl>    <dbl>     <dbl> <chr>    
##  1 A     2014-02-28  38.1     88.5 -0.0177   High_ESG 
##  2 A     2014-03-31  37.4     88.5 -0.0313   High_ESG 
##  3 A     2014-04-30  36.3     88.5  0.0537   High_ESG 
##  4 A     2014-05-30  38.2     87.6  0.0111   High_ESG 
##  5 A     2014-06-30  38.6     87.2 -0.0235   High_ESG 
##  6 A     2014-07-31  37.7     87.2  0.0191   High_ESG 
##  7 A     2014-08-29  38.5     87.4 -0.000832 High_ESG 
##  8 A     2014-09-30  38.4     86.9 -0.0298   High_ESG 
##  9 A     2014-10-31  37.3     88.2  0.0812   High_ESG 
## 10 A     2014-11-28  40.3     88.1 -0.0421   High_ESG 
## # ā€¦ with 32,630 more rows

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

data_esg2 %>%
    group_by(Date, ESG_group) %>%
    summarise(Return = mean(Return, na.rm = T)) %>%
    group_by(ESG_group) %>%
    summarise(avg_return = mean(Return, na.rm = T),  # Be careful with na.rm = T
              vol = sd(Return, na.rm = T),
              SR = avg_return/vol)
## # A tibble: 2 Ɨ 4
##   ESG_group avg_return    vol    SR
##   <chr>          <dbl>  <dbl> <dbl>
## 1 High_ESG      0.0114 0.0463 0.246
## 2 Low_ESG       0.0107 0.0476 0.225

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_esg2 %>%
    mutate(Year = year(Date)) %>%
    group_by(Year, Tick) %>%
    summarise(avg_ESG = mean(ESG_rank, na.rm = T),
              avg_return = mean(Return, na.rm = T)) %>%
    ggplot(aes(x = avg_ESG, 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_esg2){
    data_esg2 %>%
        filter(ESG_rank > level) %>%   # Filter to keep only the firms with ESG > level
        pull(Return) %>%
        mean(na.rm = T)
}

level_esg <- 10*(1:9)                                 # ESG_rank threshold
map_dbl(level_esg, ESG_impact, data_esg2 = data_esg2) # The raw result
## [1] 0.01106186 0.01111554 0.01128353 0.01122823 0.01153817 0.01141918 0.01184463
## [8] 0.01164315 0.01184526
level_esg %>%
    bind_cols(return = map_dbl(level_esg, ESG_impact, data_esg2 = data_esg2)) %>%
    ggplot(aes(x = level_esg, y = return)) + geom_col(alpha = 0.5) + 
    theme_light() + scale_y_continuous(limits=c(0.008,0.012), oob = rescale_none)

Overall, higher ESG selection seems promising indeed (but letā€™s not jump to general conclusions). Finally, letā€™s have a look at more homogeneous groups. We need different reference points (level_esg) because of the distribution of ESG_rank.

ESG_quant <- function(level, width, data_esg2){ # A new function
    data_esg2 %>%
        filter(ESG_rank > level - width/2,
               ESG_rank < level + width/2) %>%     # Filter to keep firms with ESG around level
        pull(Return) %>%
        mean(na.rm = T)
}
level_esg <- quantile(data_esg$ESG_rank, seq(0.1,0.9, by = 0.1))
level_esg %>%
    bind_cols(return = map_dbl(level_esg, ESG_quant, width = 10, data_esg2 = data_esg2)) %>%
    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.