Five measures of stock market valuation

################################################################################
###
### A dashboard for stock market valuation.   Uses the following metrics:
###
###  * Shiller PE (or CAPE)
###  * Tobin's Q
###  * AIEA (Average Investor Equity Allocation)
###  * "Buffett Ratio"
###  * Detrended log(stock price)
###
################################################################################

### You'll need to get a FRED API key from the URL below to access FRED.
### https://research.stlouisfed.org/useraccount/apikeys

api_key_fred <- "PUT_YOUR_FRED_API_KEY_HERE"
fredr_set_key(api_key_fred)

### Path to whereever you saved the Tobin's Q historical data
tobins_q_historical_csv <- "tobins_q_historical_data.csv"

################################################################################

library(fredr)
library(tidyverse)
library(cowplot)
library(lubridate)
library(tsibble)
library(broom)
library(scales)

################################################################################
### Function to add recession bars to ggplot graphs
################################################################################

geom_recession_bars <- function(date_start, date_end, fill = "darkseagreen4") {

  date_start <- as.Date(date_start, origin = "1970-01-01")
  date_end   <- as.Date(date_end,   origin = "1970-01-01")

  recessions_tibble <-
    tribble(
      ~peak,     ~trough,
      "1857-06-01", "1858-12-01",
      "1860-10-01", "1861-06-01",
      "1865-04-01", "1867-12-01",
      "1869-06-01", "1870-12-01",
      "1873-10-01", "1879-03-01",
      "1882-03-01", "1885-05-01",
      "1887-03-01", "1888-04-01",
      "1890-07-01", "1891-05-01",
      "1893-01-01", "1894-06-01",
      "1895-12-01", "1897-06-01",
      "1899-06-01", "1900-12-01",
      "1902-09-01", "1904-08-01",
      "1907-05-01", "1908-06-01",
      "1910-01-01", "1912-01-01",
      "1913-01-01", "1914-12-01",
      "1918-08-01", "1919-03-01",
      "1920-01-01", "1921-07-01",
      "1923-05-01", "1924-07-01",
      "1926-10-01", "1927-11-01",
      "1929-08-01", "1933-03-01",
      "1937-05-01", "1938-06-01",
      "1945-02-01", "1945-10-01",
      "1948-11-01", "1949-10-01",
      "1953-07-01", "1954-05-01",
      "1957-08-01", "1958-04-01",
      "1960-04-01", "1961-02-01",
      "1969-12-01", "1970-11-01",
      "1973-11-01", "1975-03-01",
      "1980-01-01", "1980-07-01",
      "1981-07-01", "1982-11-01",
      "1990-07-01", "1991-03-01",
      "2001-03-01", "2001-11-01",
      "2007-12-01", "2009-06-01",
      "2020-02-01", "2020-04-01"
    ) %>%
    mutate(peak   = as.Date(peak),
           trough = as.Date(trough))

  recessions_trim <- recessions_tibble %>%
    filter(peak   >= min(date_start) &
             trough <= max(date_end))

  if (nrow(recessions_trim) > 0) {

    recession_bars <- geom_rect(data        = recessions_trim,
                                inherit.aes = F,
                                fill        = fill,
                                alpha       = 0.25,
                                aes(xmin = as.Date(peak,   origin = "1970-01-01"),
                                    xmax = as.Date(trough, origin = "1970-01-01"),
                                    ymin = -Inf, ymax = +Inf))
  } else {
    recession_bars <- geom_blank()
  }
}

################################################################################
### Functions for reading data from .XLS/.XLSX files
################################################################################

read_excel_url <- function(url, ...) {
  tf <- tempfile(fileext = ".xls")
  curl::curl_download(url, tf)
  return(readxl::read_excel(tf, ...))
}

read_xlsx_url <- function(url, ...) {
  tf <- tempfile(fileext = ".xlsx")
  curl::curl_download(url, tf)
  return(readxl::read_excel(tf, ...))
}

################################################################################
### Pull CAPE from Shiller's spreadsheet
################################################################################

spreadsheet <-
  read_excel_url("http://www.econ.yale.edu/~shiller/data/ie_data.xls",
                 sheet = "Data", skip = 9,
                 col_names = c("Date_Foo",           "SP_Comp",
                               "Dividend",           "Earnings",
                               "CPI",                "Date_Fraction",
                               "Long_Interest_Rate", "Real_Price",
                               "Real_Dividend",      "Real_Total_Return",
                               "Real_Earnings",      "Real_TR_Scaled_Earnings",
                               "CAPE",               "Unknown1",
                               "TR_CAPE",            "Unknown2",
                               "Excess_CAPE_Yield",
                               "Monthly_Total_Bond_Returns",
                               "Real_Total_Bond_Returns",
                               "10Yr_Annualized_Stock_Real_Return",
                               "10Yr_Annualized_Bond_Real_Return",
                               "Real_10_Year_Excess_Annualized_Returns"))

data_shiller_pe <-
  spreadsheet %>%
  select(Date_Fraction, CAPE) %>%
  filter(CAPE != "NA") %>%
  mutate(CAPE = as.numeric(CAPE),
         Date = date_decimal(Date_Fraction),
         YearMonth = yearmonth(Date)) %>%
  select(YearMonth, CAPE)

################################################################################
### Compute the SP500 mean reversion
################################################################################

sp_500 <-
  spreadsheet %>%
  rename(SP_Price = Real_Price) %>%
  mutate(log_SP_Price = log(SP_Price)) %>%
  select(Date_Fraction, SP_Price, log_SP_Price) %>%
  filter(!is.na(SP_Price)) %>%
  mutate(YearMonth = yearmonth(date_decimal(Date_Fraction))) %>%
  select(YearMonth, Date_Fraction, SP_Price, log_SP_Price)

fit_mr <-
  nls(log(SP_Price) ~ a + b * Date_Fraction,
      start = c(a = -30, b = 0.0190),
      data = sp_500) %>%
  tidy()

data_sp500_mean_reversion <-
  sp_500 %>%
  mutate(
    fit_log_SP_Price = fit_mr$estimate[1] + fit_mr$estimate[2] * Date_Fraction,
    SP_Price_Trend = exp(fit_log_SP_Price),
    SP_Mean_Reversion_Delta = log_SP_Price - fit_log_SP_Price
    )

################################################################################
### Compute Aggregate Investor Equity Allocation
################################################################################

series_aiea <-
  c("NCBEILQ027S", # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
    "FBCELLQ027S", # Domestic Financial Sectors; Corporate Equities; Liability, Level
    "BCNSDODNS",   # Nonfinancial Corporate Business; Debt Securities and Loans; Liability, Level
    "CMDEBT",      # Households and Nonprofit Organizations; Debt Securities and Loans; Liability, Level

    "FGSDODNS",    # Federal Government; Debt Securities and Loans; Liability, Level
    "SLGSDODNS",   # State and Local Governments; Debt Securities and Loans; Liability, Level
    "DODFFSWCMI"   # Rest of the World; Debt Securities and Loans; Liability, Level
  )

### Use purrr to pull all of the series from FRED
data_aiea <-
  purrr::map_dfr(series_aiea, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  select(date, series_id, value) %>%
  pivot_wider(id_cols = "date",
              names_from = "series_id",
              values_from = "value") %>%
  arrange(date) %>%
  filter(date >= as.Date("1951-10-01")) %>%
  mutate(YearMonth = yearmonth(date)) %>%
  mutate(Equity = (NCBEILQ027S + FBCELLQ027S) / 1000,
         Debt   = BCNSDODNS + CMDEBT + FGSDODNS + SLGSDODNS + DODFFSWCMI,
         AIEA   = Equity / (Equity + Debt)) %>%
  select(YearMonth, AIEA)


################################################################################
### Compute the Buffett Ratio
################################################################################
series_buffett <-
  c(
    "WILL5000PRFC",
    "GDP",
    "WALCL",
    "NCBEILQ027S" # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
  )

data_buffett <-
  purrr::map_dfr(series_buffett, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  select(date, series_id, value) %>%
  pivot_wider(id_cols = "date",
              names_from = "series_id",
              values_from = "value") %>%
  arrange(date) %>%
  mutate(NCBEILQ027S = NCBEILQ027S / 1000) %>%

  # We use the alternative series whenever Wilshire 5000 isn't available from FRED
  #mutate(numerator = if_else(is.na(WILL5000PRFC), NCBEILQ027S, WILL5000PRFC)) %>%
  mutate(numerator = NCBEILQ027S) %>%

  filter(!is.na(GDP),
         !is.na(numerator)) %>%
  mutate(Buffett_Indicator = if_else(!is.na(WALCL),
                                     100 * numerator / (GDP + WALCL/1000),
                                     100 * numerator / (GDP             ))) %>%
  mutate(YearMonth = yearmonth(date)) %>%
  select(YearMonth, Buffett_Indicator)

# For data prior to 1970 (where Wilshire data is not available) we use the
# Z.1 Financial Account - Nonfinancial corporate business; corporate equities;
# liability, Level, published by the Federal Reserve, which provides a quarterly
# estimate of total market value back to 1945. In order to integrate the
# datasets, we index the Z.1 data to match up to the 1970 Wilshire starting point.
#
# Ref:  https://www.currentmarketvaluation.com/models/buffett-indicator.php

################################################################################
### Compute Tobin's Q
################################################################################

# Historical data taken from http://www.smithers.co.uk/download.php?file=918
data_tobins_q_historical <-
  read_csv(tobins_q_historical_csv) %>%
  mutate(quarter = case_when(
    quarter == "Q1" ~ 0.00,
    quarter == "Q2" ~ 0.25,
    quarter == "Q3" ~ 0.50,
    quarter == "Q4" ~ 0.75,
  )) %>%
  mutate(decimal_year = year + quarter) %>%
  mutate(date = date_decimal(decimal_year)) %>%
  mutate(YearMonth = yearmonth(date)) %>%
  rename(Tobins_Q = tobins_q) %>%
  select(YearMonth, Tobins_Q)

# Ref:  https://www.advisorperspectives.com/dshort/updates/2021/07/06/the-q-ratio-and-market-valuation-june-update

series_tobins_q <-
  c("NCBEILQ027S", # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
    "TNWMVBSNNCB"  # Nonfinancial Corporate Business; Net Worth, Level
  )

data_tobins_q <-
  purrr::map_dfr(series_tobins_q, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  select(date, series_id, value) %>%
  pivot_wider(id_cols = "date",
              names_from = "series_id",
              values_from = "value") %>%
  arrange(date) %>%
  filter(date >= as.Date("1951-10-01")) %>%
  mutate(Tobins_Q = NCBEILQ027S / (1000 * TNWMVBSNNCB)) %>%
  select(date, Tobins_Q) %>%
  mutate(YearMonth = yearmonth(date)) %>%
  select(YearMonth, Tobins_Q)

data_tobins_q <-
  data_tobins_q_historical %>%
  bind_rows(data_tobins_q) %>%
  arrange(YearMonth) %>%
  unique()

################################################################################
### Merge data and compute mean/standard deviation for each metric
################################################################################

data_valuation <-
  data_shiller_pe %>%
  left_join(data_sp500_mean_reversion, by = "YearMonth") %>%
  left_join(data_tobins_q, by = "YearMonth") %>%
  left_join(data_aiea, by = "YearMonth") %>%
  left_join(data_buffett, by = "YearMonth") %>%
  mutate(Date = as.Date(YearMonth)) %>%
  select(Date, YearMonth, CAPE, Tobins_Q, AIEA, Buffett_Indicator, SP_Mean_Reversion_Delta)

mean_shiller_pe <- data_valuation %>% pull(CAPE) %>% mean()
sd_shiller_pe   <- data_valuation %>% pull(CAPE) %>% sd()

mean_aiea       <- data_valuation %>% filter(!is.na(AIEA)) %>% pull(AIEA) %>% mean()
sd_aiea         <- data_valuation %>% filter(!is.na(AIEA)) %>% pull(AIEA) %>% sd()

mean_buffett    <- data_valuation %>% filter(!is.na(Buffett_Indicator)) %>% pull(Buffett_Indicator) %>% mean()
sd_buffett      <- data_valuation %>% filter(!is.na(Buffett_Indicator)) %>% pull(Buffett_Indicator) %>% sd()

mean_tobins_q   <- data_valuation %>% filter(!is.na(Tobins_Q)) %>% pull(Tobins_Q) %>% mean()
sd_tobins_q     <- data_valuation %>% filter(!is.na(Tobins_Q)) %>% pull(Tobins_Q) %>% sd()

mean_sp500_mr   <- data_valuation %>% filter(!is.na(SP_Mean_Reversion_Delta)) %>% pull(SP_Mean_Reversion_Delta) %>% mean()
sd_sp500_mr     <- data_valuation %>% filter(!is.na(SP_Mean_Reversion_Delta)) %>% pull(SP_Mean_Reversion_Delta) %>% sd()

################################################################################
### Assemble graphs...
################################################################################

g_shiller_pe_z <-
  ggplot(data = data_valuation) +
  theme_bw() +
  geom_line(aes(x = Date, y = CAPE)) +

  geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +

  geom_hline(yintercept = mean_shiller_pe + 3 * sd_shiller_pe, linetype = "dotted", color = "darkred") +
  geom_hline(yintercept = mean_shiller_pe + 2 * sd_shiller_pe, linetype = "dotted", color = "red") +
  geom_hline(yintercept = mean_shiller_pe + 1 * sd_shiller_pe, linetype = "dotted", color = "firebrick1") +
  geom_hline(yintercept = mean_shiller_pe + 0 * sd_shiller_pe, linetype = "dashed", color = "black") +
  geom_hline(yintercept = mean_shiller_pe - 1 * sd_shiller_pe, linetype = "dotted", color = "blue") +
  geom_hline(yintercept = mean_shiller_pe - 2 * sd_shiller_pe, linetype = "dotted", color = "darkblue") +

  labs(x = "",
       y = "",
       #y = "Shiller PE",
       title = "Shiller PE / CAPE") +
  scale_x_date(breaks = pretty_breaks(9)) +
  scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
                     sec.axis = sec_axis(~ (. - mean_shiller_pe) / sd_shiller_pe,
                                         name = "Z-Score",
                                         breaks = pretty_breaks(6)))
print(g_shiller_pe_z)


g_sp500_mr <-
  ggplot(data = data_valuation) +
  theme_bw() +
  geom_line(aes(x = Date, y = SP_Mean_Reversion_Delta)) +

  geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +

  geom_hline(yintercept = mean_sp500_mr + 3 * sd_sp500_mr, linetype = "dotted", color = "darkred") +
  geom_hline(yintercept = mean_sp500_mr + 2 * sd_sp500_mr, linetype = "dotted", color = "red") +
  geom_hline(yintercept = mean_sp500_mr + 1 * sd_sp500_mr, linetype = "dotted", color = "firebrick1") +
  geom_hline(yintercept = mean_sp500_mr + 0 * sd_sp500_mr, linetype = "dashed", color = "black") +
  geom_hline(yintercept = mean_sp500_mr - 1 * sd_sp500_mr, linetype = "dotted", color = "blue") +
  geom_hline(yintercept = mean_sp500_mr - 2 * sd_sp500_mr, linetype = "dotted", color = "darkblue") +

  labs(x = "",
       y = "",
       title = "Detrended log(Real S&P 500 Price)") +
  scale_x_date(breaks = pretty_breaks(9)) +
  scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
    sec.axis = sec_axis(~ (. - mean_sp500_mr) / sd_sp500_mr,
                        name = "Z-Score",
                        breaks = pretty_breaks(6)))
print(g_sp500_mr)


g_tobins_q_z <-
  ggplot(data = data_valuation %>% filter(!is.na(Tobins_Q))) +
  theme_bw() +
  geom_line(aes(x = Date, y = Tobins_Q)) +

  geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +

  geom_hline(yintercept = mean_tobins_q + 3 * sd_tobins_q, linetype = "dotted", color = "darkred") +
  geom_hline(yintercept = mean_tobins_q + 2 * sd_tobins_q, linetype = "dotted", color = "red") +
  geom_hline(yintercept = mean_tobins_q + 1 * sd_tobins_q, linetype = "dotted", color = "firebrick1") +
  geom_hline(yintercept = mean_tobins_q + 0 * sd_tobins_q, linetype = "dashed", color = "black") +
  geom_hline(yintercept = mean_tobins_q - 1 * sd_tobins_q, linetype = "dotted", color = "blue") +
  geom_hline(yintercept = mean_tobins_q - 2 * sd_tobins_q, linetype = "dotted", color = "darkblue") +

  labs(x = "",
       y = "",
       #y = "Tobin's Q",
       title = "Tobin's Q") +
  scale_x_date(breaks = pretty_breaks(9)) +
  scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
    sec.axis = sec_axis(~ (. - mean_tobins_q) / sd_tobins_q,
                        name = "Z-Score",
                        breaks = pretty_breaks(6)))
print(g_tobins_q_z)


g_aiea_z <-
  ggplot(data = data_valuation %>% filter(!is.na(AIEA))) +
  theme_bw() +
  geom_line(aes(x = Date, y = AIEA)) +

  geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +

  geom_hline(yintercept = mean_aiea + 2 * sd_aiea, linetype = "dotted", color = "red") +
  geom_hline(yintercept = mean_aiea + 1 * sd_aiea, linetype = "dotted", color = "firebrick1") +
  geom_hline(yintercept = mean_aiea + 0 * sd_aiea, linetype = "dashed", color = "black") +
  geom_hline(yintercept = mean_aiea - 1 * sd_aiea, linetype = "dotted", color = "blue") +
  geom_hline(yintercept = mean_aiea - 2 * sd_aiea, linetype = "dotted", color = "darkblue") +

  labs(x = "",
       y = "",
       #y = "AIEA",
       title = "Aggregate Investor Equity Allocation") +
  scale_x_date(breaks = pretty_breaks(9)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
    sec.axis = sec_axis(~ (. - mean_aiea) / sd_aiea,
                        name = "Z-Score",
                        breaks = pretty_breaks(6)))
print(g_aiea_z)

g_buffett_z <-
  ggplot(data = data_valuation %>% filter(!is.na(Buffett_Indicator))) +
  theme_bw() +
  geom_line(aes(x = Date, y = Buffett_Indicator)) +

  geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +

  geom_hline(yintercept = mean_buffett + 3 * sd_buffett, linetype = "dotted", color = "darkred") +
  geom_hline(yintercept = mean_buffett + 2 * sd_buffett, linetype = "dotted", color = "red") +
  geom_hline(yintercept = mean_buffett + 1 * sd_buffett, linetype = "dotted", color = "firebrick1") +
  geom_hline(yintercept = mean_buffett + 0 * sd_buffett, linetype = "dashed", color = "black") +
  geom_hline(yintercept = mean_buffett - 1 * sd_buffett, linetype = "dotted", color = "blue") +
  geom_hline(yintercept = mean_buffett - 2 * sd_buffett, linetype = "dotted", color = "darkblue") +

  labs(x = "",
       y = "",
       title = "Buffett Indicator (Total Stock Market Cap / [GDP + Total Fed Assets])") +
  scale_x_date(breaks = pretty_breaks(9)) +
  scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
                     sec.axis = sec_axis(~ (. - mean_buffett) / sd_buffett,
                                         name = "Z-Score",
                                         breaks = pretty_breaks(6)))
print(g_buffett_z)

print(
  plot_grid(
    g_shiller_pe_z,
    g_tobins_q_z,
    g_aiea_z,
    g_buffett_z,
    g_sp500_mr,
    nrow = 5, ncol = 1, align = "hv"))