Welcome to another installment of Reproducible Finance. Today’s post will be topical as we look at the historical behavior of the stock market after days of extreme returns and it will also explore one of my favorite coding themes of 2020 – the power of RMarkdown as an R/Python collaboration tool.

This post originated when Rishi Singh, the founder of tiingo and one of the nicest people I have encountered in this crazy world, sent over a note about recent market volatility along with some Python code for analyzing that volatility. We thought it would be a nice project to post that Python code along with the equivalent R code for reproducing the same results. For me, it’s a great opportunity to use RMarkdown's R and Python interoperability superpowers, fueled by the reticulate package. If you are an R coder and someone sends you Python code as part of a project, RMarkdown + reticulate makes it quite smooth to incorporate that Python code into your work. It was interesting to learn how a very experienced Python coder might tackle a problem and then think about how to tackle that problem with R. Unsurprisingly, I couldn’t resist adding a few elements of data visualization.

Before we get started, if you’re unfamiliar with using R and Python chunks throughout an RMarkdown file, have a quick look at the reticulate documentation here.

Let’s get to it. Since we’ll be working with R and Python, we start with our usual R setup code chunk to load R packages, but we’ll also load the reticulate package and source a Python script. Here’s what that looks like.

library(tidyverse)
library(tidyquant)
library(riingo)
library(timetk)
library(plotly)
library(roll)
library(slider)
library(reticulate)

riingo_set_token("your tiingo token here")

# Python file that holds my tiingo token
reticulate::source_python("credentials.py")

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = NA)

Note that I set my tiingo token twice: first using riingo_set_token() so I can use the riingo package in R chunks and then by sourcing the credentials.py file, where I have put tiingoToken = 'my token'. Now I can use the tiingoToken variable in my Python chunks. This is necessary because we will use both R and Python to pull in data from Tiingo.

Next we will use a Python chunk to load the necessary Python libraries. If you haven’t installed these yet, you can open the RStudio terminal and run pip install. Since we’ll be interspersing R and Python code chunks throughout, I will add a # Python Chunk to each Python chunk and, um, # R Chunk to each R chunk.

# Python chunk
import pandas as pd
import numpy as np
import tiingo

Let’s get to the substance. The goal today is look back at the last 43 years of S&P 500 price history and analyze how the market has performed following a day that sees an extreme return. We will also take care with how we define an extreme return, using rolling volatility to normalize percentage moves.

We will use the mutual fund VFINX as a tradeable proxy for the S&P 500 because it has a much longer history than other funds like SPY or VOO.

Let’s start by passing a URL string from tiingo to the pandas function read_csv, along with our tiingoToken.

# Python chunk 
pricesDF = pd.read_csv("https://api.tiingo.com/tiingo/daily/vfinx/prices?startDate=1976-1-1&format=csv&token=" + tiingoToken)

We just created a Python object called pricesDF. We can look at that object in an R chunk by calling py$pricesDF.

# R chunk
py$pricesDF %>% 
  head()

We just created a Python object called pricesDF. Let’s reformat the date column becomes the index, in date time format.

# Python chunk 
pricesDF = pricesDF.set_index(['date'])
pricesDF.index = pd.DatetimeIndex(pricesDF.index)

Heading back to R for viewing, we see that the date column is no longer a column – it is the index of the data frame and in pandas the index is more like a label than a new column. In fact, here’s what happens when call the row names of this data frame.

# R chunk
py$pricesDF %>% 
  head() %>% 
  rownames()

We now have our prices, indexed by date. Let’s convert adjusted closing prices to log returns and save the results in a new column called returns. Note the use of the shift(1) operator here. That is analogous to the lag(..., 1) function in dplyr.

# Python chunk
pricesDF['returns'] = np.log(pricesDF['adjClose']/pricesDF['adjClose'].shift(1))

Next, we want to calculate the 3-month rolling standard deviation of these daily log returns, and then divide daily returns by the previous rolling 3-month volatility in order to prevent look-ahead error. We can think of this as normalizing today’s return by the previous 3-months’ rolling volatility and will label it as stdDevMove.

# Python chunk
pricesDF['rollingVol'] = pricesDF['returns'].rolling(63).std()
pricesDF['stdDevMove'] = pricesDF['returns'] / pricesDF['rollingVol'].shift(1)

Finally, we eventually want to calculate how the market has performed on the day following a large negative move. To prepare for that, let’s create a column of next day returns using shift(-1).

# Python chunk
pricesDF['nextDayReturns'] = pricesDF.returns.shift(-1)

Now, we can filter by the size of the stdDevMove column and the returns column, to isolate days where the standard deviation move was at least 3 and the returns was less than -3%. We use mean() to find the mean next day return following such large events.

# Python chunk
nextDayPerformanceSeries = pricesDF.loc[(pricesDF['stdDevMove'] < -3) & (pricesDF['returns'] < -.03), ['nextDayReturns']].mean()

Finally, let’s loop through and see how the mean next day return changes as we filter on different extreme negative returns or we can call drop tolerances. We will label the drop tolerance as i, set it at -.03 and then run a while loop that decrements down i by .0025 at each pass. In this way we can look at the mean next return following different levels of negative returns.

# Python chunk
i = -.03
while i >= -.0525:
    nextDayPerformanceSeries = pricesDF.loc[(pricesDF['stdDevMove'] < -3) & (pricesDF['returns'] < i), ['nextDayReturns']]
    print(str(round(i, 5)) + ': ' + str(round(nextDayPerformanceSeries['nextDayReturns'].mean(), 6)))

    i -= .0025

It appears that as the size of the drop gets larger and more negative, the mean bounce back tends to get larger.

Let’s reproduce these results in R.

First, we import prices using the riingo_prices() function from the riingo.

# R chunk
sp_500_prices <- 
  "VFINX" %>% 
  riingo_prices(start_date = "1976-01-01", end_date = today())

We can use mutate() to add a column of daily returns, rolling volatility, standard deviation move and next day returns.

# R chunk
sp_500_returns <- 
  sp_500_prices %>%
  select(date, adjClose) %>% 
  mutate(daily_returns_log = log(adjClose/lag(adjClose)),
         rolling_vol = roll_sd(as.matrix(daily_returns_log), 63),
         sd_move = daily_returns_log/lag(rolling_vol),
         next_day_returns = lead(daily_returns_log))

Now let’s filter() on an sd_move greater than 3 and daily_returns_log less than a drop tolerance of -.03.

# R chunk

sp_500_returns %>% 
  na.omit() %>% 
  filter(sd_move < -3 & daily_returns_log < -.03) %>% 
  select(date, daily_returns_log, sd_move, next_day_returns) %>% 
  summarise(mean_return = mean(next_day_returns)) %>% 
  add_column(drop_tolerance = scales::percent(.03), .before = 1)
# A tibble: 1 x 2
  drop_tolerance mean_return
  <chr>                <dbl>
1 3%                 0.00625

We used a while() loop to iterate across different drop tolerances in Python, let’s see how to implement that using map_dfr() from the purrr package.

First, we will define a sequence of drop tolerances using the seq() function.

# R chunk
drop_tolerance <- seq(.03, .05, .0025)

drop_tolerance
[1] 0.0300 0.0325 0.0350 0.0375 0.0400 0.0425 0.0450 0.0475 0.0500

Next, we will create a function called outlier_mov_fun that takes a data frame of returns, filters on a drop tolerance and gives us the mean return following large negative moves.

# R chunk
outlier_mov_fun <- function(drop_tolerance, returns) {
returns %>%  
  na.omit() %>% 
  filter(sd_move < -3 & daily_returns_log < -drop_tolerance) %>% 
  select(date, daily_returns_log, sd_move, next_day_returns) %>% 
  summarise(mean_return = mean(next_day_returns) %>% round(6)) %>%
  add_column(drop_tolerance = scales::percent(drop_tolerance), .before = 1) %>% 
  add_column(drop_tolerance_raw = drop_tolerance, .before = 1)
}

Notice how that function takes two arguments: a drop tolerance and data frame of returns.

Next, we pass our sequence of drop tolerances, stored in a variable called drop_tolerance to map_dfr(), along with our function and our sp_500_returns object. map_dfr will iterate through our sequence of drops and apply our function to each one.

# R chunk
map_dfr(drop_tolerance, outlier_mov_fun, sp_500_returns) %>% 
  select(-drop_tolerance_raw)
# A tibble: 9 x 2
  drop_tolerance mean_return
  <chr>                <dbl>
1 3%                 0.00625
2 3%                 0.00700
3 3%                 0.00967
4 4%                 0.0109 
5 4%                 0.0122 
6 4%                 0.0132 
7 4%                 0.0149 
8 5%                 0.0149 
9 5%                 0.0162 

Have a quick glance up that the results of our Python while() and we should see that the results are consistent.

Alright, let’s have some fun and get to visualizing these results with ggplot and plotly.

# R chunk
(
sp_500_returns %>% 
map_dfr(drop_tolerance, outlier_mov_fun, .) %>% 
  ggplot(aes(x = drop_tolerance_raw, y = mean_return, text = str_glue("drop tolerance: {drop_tolerance}
                                                                      mean next day return: {mean_return * 100}%"))) +
  geom_point(color = "cornflowerblue") +
  labs(title = "Mean Return after Large Daily Drop", y = "mean return", x = "daily drop") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::percent) + 
  theme_minimal()
) %>% ggplotly(tooltip = "text")

Here’s what happens when we expand the upper bound to a drop tolerance of -2% and make our intervals smaller, moving from .25% increments to .125% increments.

# R chunk
drop_tolerance_2 <- seq(.02, .05, .00125)

(
sp_500_returns %>% 
map_dfr(drop_tolerance_2, outlier_mov_fun, .) %>% 
  ggplot(aes(x = drop_tolerance_raw, y = mean_return, text = str_glue("drop tolerance: {drop_tolerance}
                                                                      mean next day return: {mean_return * 100}%"))) +
  geom_point(color = "cornflowerblue") +
  labs(title = "Mean Return after Large Daily Drop", y = "mean return", x = "daily drop") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::percent) + 
  theme_minimal()
) %>% ggplotly(tooltip = "text")

Check out what happens when we expand the lower bound, to a -6% drop tolerance.

# R chunk
drop_tolerance_3 <- seq(.02, .06, .00125)

(
sp_500_returns %>% 
map_dfr(drop_tolerance_3, outlier_mov_fun, .) %>% 
  ggplot(aes(x = drop_tolerance_raw, y = mean_return, text = str_glue("drop tolerance: {drop_tolerance}
                                                                      mean next day return: {mean_return * 100}%"))) +
  geom_point(color = "cornflowerblue") +
  labs(title = "Mean Return after Large Daily Drop", y = "mean return", x = "daily drop") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::percent) + 
  theme_minimal()
) %>% ggplotly(tooltip = "text")

I did not expect that gap upward when the daily drop passes 5.25%.

A quick addendum that if I had gotten my act together and finished this 4 days ago I would not have included, but I’m curious how this last week has compared with other weeks in terms of volatility. I have in mind to visualize weekly return dispersion and that seemed a mighty tall task, until the brand new slider package came to the rescue! slider has a function called slide_period() that, among other things, allows us to break up time series according to different periodicities.

To break up our returns by week, we call slide_period_dfr(., .$date, "week", ~ .x, .origin = first_monday_december, .names_to = "week"), where first_monday_december is a date that falls on a Monday. We could use our eyeballs to check a calendar and find a date that’s a Monday or we could use some good ol’ code. Let’s assume we want to find the first Monday in December of 2016.

We first filter our data with filter(between(date, as_date("2016-12-01"), as_date("2016-12-31"))). Then create a column of weekday names with wday(date, label = TRUE, abbr = FALSE) and filter to our first value of “Monday”.

# R Chunk
first_monday_december <- 
  sp_500_returns %>%
  mutate(date = ymd(date)) %>% 
  filter(between(date, as_date("2016-12-01"), as_date("2016-12-31"))) %>% 
  mutate(day_week = wday(date, label = TRUE, abbr = FALSE)) %>% 
  filter(day_week == "Monday") %>% 
  slice(1) %>% 
  pull(date)

Now we run our slide_period_dfr() code and it will start on the first Monday in December of 2016, and break our returns into weeks. Since we set .names_to = "week", the function will create a new column called week and give a unique number to each of our weeks.

# R chunk

sp_500_returns %>%
  select(date, daily_returns_log) %>%
  filter(date >= first_monday_december) %>%
  slide_period_dfr(.,
                   .$date,
                   "week",
                   ~ .x,
                   .origin = first_monday_december,
                   .names_to = "week") %>% 
  head(10)
# A tibble: 10 x 3
    week date                daily_returns_log
   <int> <dttm>                          <dbl>
 1     1 2016-12-05 00:00:00           0.00589
 2     1 2016-12-06 00:00:00           0.00342
 3     1 2016-12-07 00:00:00           0.0133 
 4     1 2016-12-08 00:00:00           0.00226
 5     1 2016-12-09 00:00:00           0.00589
 6     2 2016-12-12 00:00:00          -0.00105
 7     2 2016-12-13 00:00:00           0.00667
 8     2 2016-12-14 00:00:00          -0.00810
 9     2 2016-12-15 00:00:00           0.00392
10     2 2016-12-16 00:00:00          -0.00172

From here, we can group_by that week column and treat each week as a discrete time period. Let’s use ggplotly to plot each week on the x-axis and the daily returns of each week on the y-axis, so that the vertical dispersion shows us the dispersion of weekly returns. Hover on the point to see the exact date of the return.

# R chunk
(
sp_500_returns %>%
  select(date, daily_returns_log) %>%
  filter(date >= first_monday_december) %>%
  slide_period_dfr(.,
                   .$date,
                   "week",
                   ~ .x,
                   .origin = first_monday_december,
                   .names_to = "week") %>%
  group_by(week) %>%
  mutate(start_week = ymd(min(date))) %>%
  ggplot(aes(x = start_week, y = daily_returns_log, text = str_glue("date: {date}"))) +
  geom_point(color = "cornflowerblue", alpha = .5) +
  scale_y_continuous(labels = scales::percent,
                     breaks = scales::pretty_breaks(n = 8)) +
  scale_x_date(breaks = scales::pretty_breaks(n = 10)) +
  labs(y = "", x = "", title = "Weekly Daily Returns") +
  theme_minimal()
) %>% ggplotly(tooltip = "text")

We can also plot the standard deviation of returns for each week.

# R chunk

(
sp_500_returns %>%
  select(date, daily_returns_log) %>%
  filter(date >= first_monday_december) %>%
  slide_period_dfr(.,
                   .$date,
                   "week",
                   ~ .x,
                   .origin = first_monday_december,
                   .names_to = "week") %>%
  group_by(week) %>%
  summarise(first_of_week = first(date),
            sd = sd(daily_returns_log)) %>%
  ggplot(aes(x = first_of_week, y = sd, text = str_glue("week: {first_of_week}"))) +
  geom_point(aes(color = sd)) +
  labs(x = "", title = "Weekly Standard Dev of Returns", y = "") +
  theme_minimal()
) %>% ggplotly(tooltip = "text")

That’s all for today! Thanks for reading and stay safe out there.