Advertisement
MetricT

Five measures of stock market valuation

Aug 11th, 2021 (edited)
2,800
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 18.78 KB
  1. ################################################################################
  2. ###
  3. ### A dashboard for stock market valuation.   Uses the following metrics:
  4. ###
  5. ###  * Shiller PE (or CAPE)
  6. ###  * Tobin's Q
  7. ###  * AIEA (Average Investor Equity Allocation)
  8. ###  * "Buffett Ratio"
  9. ###  * Detrended log(stock price)
  10. ###
  11. ################################################################################
  12.  
  13. ### You'll need to get a FRED API key from the URL below to access FRED.
  14. ### https://research.stlouisfed.org/useraccount/apikeys
  15.  
  16. api_key_fred <- "PUT_YOUR_FRED_API_KEY_HERE"
  17. fredr_set_key(api_key_fred)
  18.  
  19. ### Path to whereever you saved the Tobin's Q historical data
  20. tobins_q_historical_csv <- "tobins_q_historical_data.csv"
  21.  
  22. ################################################################################
  23.  
  24. library(fredr)
  25. library(tidyverse)
  26. library(cowplot)
  27. library(lubridate)
  28. library(tsibble)
  29. library(broom)
  30. library(scales)
  31.  
  32. ################################################################################
  33. ### Function to add recession bars to ggplot graphs
  34. ################################################################################
  35.  
  36. geom_recession_bars <- function(date_start, date_end, fill = "darkseagreen4") {
  37.  
  38.   date_start <- as.Date(date_start, origin = "1970-01-01")
  39.   date_end   <- as.Date(date_end,   origin = "1970-01-01")
  40.  
  41.   recessions_tibble <-
  42.     tribble(
  43.       ~peak,     ~trough,
  44.       "1857-06-01", "1858-12-01",
  45.       "1860-10-01", "1861-06-01",
  46.       "1865-04-01", "1867-12-01",
  47.       "1869-06-01", "1870-12-01",
  48.       "1873-10-01", "1879-03-01",
  49.       "1882-03-01", "1885-05-01",
  50.       "1887-03-01", "1888-04-01",
  51.       "1890-07-01", "1891-05-01",
  52.       "1893-01-01", "1894-06-01",
  53.       "1895-12-01", "1897-06-01",
  54.       "1899-06-01", "1900-12-01",
  55.       "1902-09-01", "1904-08-01",
  56.       "1907-05-01", "1908-06-01",
  57.       "1910-01-01", "1912-01-01",
  58.       "1913-01-01", "1914-12-01",
  59.       "1918-08-01", "1919-03-01",
  60.       "1920-01-01", "1921-07-01",
  61.       "1923-05-01", "1924-07-01",
  62.       "1926-10-01", "1927-11-01",
  63.       "1929-08-01", "1933-03-01",
  64.       "1937-05-01", "1938-06-01",
  65.       "1945-02-01", "1945-10-01",
  66.       "1948-11-01", "1949-10-01",
  67.       "1953-07-01", "1954-05-01",
  68.       "1957-08-01", "1958-04-01",
  69.       "1960-04-01", "1961-02-01",
  70.       "1969-12-01", "1970-11-01",
  71.       "1973-11-01", "1975-03-01",
  72.       "1980-01-01", "1980-07-01",
  73.       "1981-07-01", "1982-11-01",
  74.       "1990-07-01", "1991-03-01",
  75.       "2001-03-01", "2001-11-01",
  76.       "2007-12-01", "2009-06-01",
  77.       "2020-02-01", "2020-04-01"
  78.     ) %>%
  79.     mutate(peak   = as.Date(peak),
  80.            trough = as.Date(trough))
  81.  
  82.   recessions_trim <- recessions_tibble %>%
  83.     filter(peak   >= min(date_start) &
  84.              trough <= max(date_end))
  85.  
  86.   if (nrow(recessions_trim) > 0) {
  87.  
  88.     recession_bars <- geom_rect(data        = recessions_trim,
  89.                                 inherit.aes = F,
  90.                                 fill        = fill,
  91.                                 alpha       = 0.25,
  92.                                 aes(xmin = as.Date(peak,   origin = "1970-01-01"),
  93.                                     xmax = as.Date(trough, origin = "1970-01-01"),
  94.                                     ymin = -Inf, ymax = +Inf))
  95.   } else {
  96.     recession_bars <- geom_blank()
  97.   }
  98. }
  99.  
  100. ################################################################################
  101. ### Functions for reading data from .XLS/.XLSX files
  102. ################################################################################
  103.  
  104. read_excel_url <- function(url, ...) {
  105.   tf <- tempfile(fileext = ".xls")
  106.   curl::curl_download(url, tf)
  107.   return(readxl::read_excel(tf, ...))
  108. }
  109.  
  110. read_xlsx_url <- function(url, ...) {
  111.   tf <- tempfile(fileext = ".xlsx")
  112.   curl::curl_download(url, tf)
  113.   return(readxl::read_excel(tf, ...))
  114. }
  115.  
  116. ################################################################################
  117. ### Pull CAPE from Shiller's spreadsheet
  118. ################################################################################
  119.  
  120. spreadsheet <-
  121.   read_excel_url("http://www.econ.yale.edu/~shiller/data/ie_data.xls",
  122.                  sheet = "Data", skip = 9,
  123.                  col_names = c("Date_Foo",           "SP_Comp",
  124.                                "Dividend",           "Earnings",
  125.                                "CPI",                "Date_Fraction",
  126.                                "Long_Interest_Rate", "Real_Price",
  127.                                "Real_Dividend",      "Real_Total_Return",
  128.                                "Real_Earnings",      "Real_TR_Scaled_Earnings",
  129.                                "CAPE",               "Unknown1",
  130.                                "TR_CAPE",            "Unknown2",
  131.                                "Excess_CAPE_Yield",
  132.                                "Monthly_Total_Bond_Returns",
  133.                                "Real_Total_Bond_Returns",
  134.                                "10Yr_Annualized_Stock_Real_Return",
  135.                                "10Yr_Annualized_Bond_Real_Return",
  136.                                "Real_10_Year_Excess_Annualized_Returns"))
  137.  
  138. data_shiller_pe <-
  139.   spreadsheet %>%
  140.   select(Date_Fraction, CAPE) %>%
  141.   filter(CAPE != "NA") %>%
  142.   mutate(CAPE = as.numeric(CAPE),
  143.          Date = date_decimal(Date_Fraction),
  144.          YearMonth = yearmonth(Date)) %>%
  145.   select(YearMonth, CAPE)
  146.  
  147. ################################################################################
  148. ### Compute the SP500 mean reversion
  149. ################################################################################
  150.  
  151. sp_500 <-
  152.   spreadsheet %>%
  153.   rename(SP_Price = Real_Price) %>%
  154.   mutate(log_SP_Price = log(SP_Price)) %>%
  155.   select(Date_Fraction, SP_Price, log_SP_Price) %>%
  156.   filter(!is.na(SP_Price)) %>%
  157.   mutate(YearMonth = yearmonth(date_decimal(Date_Fraction))) %>%
  158.   select(YearMonth, Date_Fraction, SP_Price, log_SP_Price)
  159.  
  160. fit_mr <-
  161.   nls(log(SP_Price) ~ a + b * Date_Fraction,
  162.       start = c(a = -30, b = 0.0190),
  163.       data = sp_500) %>%
  164.   tidy()
  165.  
  166. data_sp500_mean_reversion <-
  167.   sp_500 %>%
  168.   mutate(
  169.     fit_log_SP_Price = fit_mr$estimate[1] + fit_mr$estimate[2] * Date_Fraction,
  170.     SP_Price_Trend = exp(fit_log_SP_Price),
  171.     SP_Mean_Reversion_Delta = log_SP_Price - fit_log_SP_Price
  172.     )
  173.  
  174. ################################################################################
  175. ### Compute Aggregate Investor Equity Allocation
  176. ################################################################################
  177.  
  178. series_aiea <-
  179.   c("NCBEILQ027S", # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
  180.     "FBCELLQ027S", # Domestic Financial Sectors; Corporate Equities; Liability, Level
  181.     "BCNSDODNS",   # Nonfinancial Corporate Business; Debt Securities and Loans; Liability, Level
  182.     "CMDEBT",      # Households and Nonprofit Organizations; Debt Securities and Loans; Liability, Level
  183.  
  184.     "FGSDODNS",    # Federal Government; Debt Securities and Loans; Liability, Level
  185.     "SLGSDODNS",   # State and Local Governments; Debt Securities and Loans; Liability, Level
  186.     "DODFFSWCMI"   # Rest of the World; Debt Securities and Loans; Liability, Level
  187.   )
  188.  
  189. ### Use purrr to pull all of the series from FRED
  190. data_aiea <-
  191.   purrr::map_dfr(series_aiea, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  192.   select(date, series_id, value) %>%
  193.   pivot_wider(id_cols = "date",
  194.               names_from = "series_id",
  195.               values_from = "value") %>%
  196.   arrange(date) %>%
  197.   filter(date >= as.Date("1951-10-01")) %>%
  198.   mutate(YearMonth = yearmonth(date)) %>%
  199.   mutate(Equity = (NCBEILQ027S + FBCELLQ027S) / 1000,
  200.          Debt   = BCNSDODNS + CMDEBT + FGSDODNS + SLGSDODNS + DODFFSWCMI,
  201.          AIEA   = Equity / (Equity + Debt)) %>%
  202.   select(YearMonth, AIEA)
  203.  
  204.  
  205. ################################################################################
  206. ### Compute the Buffett Ratio
  207. ################################################################################
  208. series_buffett <-
  209.   c(
  210.     "WILL5000PRFC",
  211.     "GDP",
  212.     "WALCL",
  213.     "NCBEILQ027S" # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
  214.   )
  215.  
  216. data_buffett <-
  217.   purrr::map_dfr(series_buffett, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  218.   select(date, series_id, value) %>%
  219.   pivot_wider(id_cols = "date",
  220.               names_from = "series_id",
  221.               values_from = "value") %>%
  222.   arrange(date) %>%
  223.   mutate(NCBEILQ027S = NCBEILQ027S / 1000) %>%
  224.  
  225.   # We use the alternative series whenever Wilshire 5000 isn't available from FRED
  226.   #mutate(numerator = if_else(is.na(WILL5000PRFC), NCBEILQ027S, WILL5000PRFC)) %>%
  227.   mutate(numerator = NCBEILQ027S) %>%
  228.  
  229.   filter(!is.na(GDP),
  230.          !is.na(numerator)) %>%
  231.   mutate(Buffett_Indicator = if_else(!is.na(WALCL),
  232.                                      100 * numerator / (GDP + WALCL/1000),
  233.                                      100 * numerator / (GDP             ))) %>%
  234.   mutate(YearMonth = yearmonth(date)) %>%
  235.   select(YearMonth, Buffett_Indicator)
  236.  
  237. # For data prior to 1970 (where Wilshire data is not available) we use the
  238. # Z.1 Financial Account - Nonfinancial corporate business; corporate equities;
  239. # liability, Level, published by the Federal Reserve, which provides a quarterly
  240. # estimate of total market value back to 1945. In order to integrate the
  241. # datasets, we index the Z.1 data to match up to the 1970 Wilshire starting point.
  242. #
  243. # Ref:  https://www.currentmarketvaluation.com/models/buffett-indicator.php
  244.  
  245. ################################################################################
  246. ### Compute Tobin's Q
  247. ################################################################################
  248.  
  249. # Historical data taken from http://www.smithers.co.uk/download.php?file=918
  250. data_tobins_q_historical <-
  251.   read_csv(tobins_q_historical_csv) %>%
  252.   mutate(quarter = case_when(
  253.     quarter == "Q1" ~ 0.00,
  254.     quarter == "Q2" ~ 0.25,
  255.     quarter == "Q3" ~ 0.50,
  256.     quarter == "Q4" ~ 0.75,
  257.   )) %>%
  258.   mutate(decimal_year = year + quarter) %>%
  259.   mutate(date = date_decimal(decimal_year)) %>%
  260.   mutate(YearMonth = yearmonth(date)) %>%
  261.   rename(Tobins_Q = tobins_q) %>%
  262.   select(YearMonth, Tobins_Q)
  263.  
  264. # Ref:  https://www.advisorperspectives.com/dshort/updates/2021/07/06/the-q-ratio-and-market-valuation-june-update
  265.  
  266. series_tobins_q <-
  267.   c("NCBEILQ027S", # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
  268.     "TNWMVBSNNCB"  # Nonfinancial Corporate Business; Net Worth, Level
  269.   )
  270.  
  271. data_tobins_q <-
  272.   purrr::map_dfr(series_tobins_q, .f = ~ fredr(series_id = .x, frequency = "q")) %>%
  273.   select(date, series_id, value) %>%
  274.   pivot_wider(id_cols = "date",
  275.               names_from = "series_id",
  276.               values_from = "value") %>%
  277.   arrange(date) %>%
  278.   filter(date >= as.Date("1951-10-01")) %>%
  279.   mutate(Tobins_Q = NCBEILQ027S / (1000 * TNWMVBSNNCB)) %>%
  280.   select(date, Tobins_Q) %>%
  281.   mutate(YearMonth = yearmonth(date)) %>%
  282.   select(YearMonth, Tobins_Q)
  283.  
  284. data_tobins_q <-
  285.   data_tobins_q_historical %>%
  286.   bind_rows(data_tobins_q) %>%
  287.   arrange(YearMonth) %>%
  288.   unique()
  289.  
  290. ################################################################################
  291. ### Merge data and compute mean/standard deviation for each metric
  292. ################################################################################
  293.  
  294. data_valuation <-
  295.   data_shiller_pe %>%
  296.   left_join(data_sp500_mean_reversion, by = "YearMonth") %>%
  297.   left_join(data_tobins_q, by = "YearMonth") %>%
  298.   left_join(data_aiea, by = "YearMonth") %>%
  299.   left_join(data_buffett, by = "YearMonth") %>%
  300.   mutate(Date = as.Date(YearMonth)) %>%
  301.   select(Date, YearMonth, CAPE, Tobins_Q, AIEA, Buffett_Indicator, SP_Mean_Reversion_Delta)
  302.  
  303. mean_shiller_pe <- data_valuation %>% pull(CAPE) %>% mean()
  304. sd_shiller_pe   <- data_valuation %>% pull(CAPE) %>% sd()
  305.  
  306. mean_aiea       <- data_valuation %>% filter(!is.na(AIEA)) %>% pull(AIEA) %>% mean()
  307. sd_aiea         <- data_valuation %>% filter(!is.na(AIEA)) %>% pull(AIEA) %>% sd()
  308.  
  309. mean_buffett    <- data_valuation %>% filter(!is.na(Buffett_Indicator)) %>% pull(Buffett_Indicator) %>% mean()
  310. sd_buffett      <- data_valuation %>% filter(!is.na(Buffett_Indicator)) %>% pull(Buffett_Indicator) %>% sd()
  311.  
  312. mean_tobins_q   <- data_valuation %>% filter(!is.na(Tobins_Q)) %>% pull(Tobins_Q) %>% mean()
  313. sd_tobins_q     <- data_valuation %>% filter(!is.na(Tobins_Q)) %>% pull(Tobins_Q) %>% sd()
  314.  
  315. mean_sp500_mr   <- data_valuation %>% filter(!is.na(SP_Mean_Reversion_Delta)) %>% pull(SP_Mean_Reversion_Delta) %>% mean()
  316. sd_sp500_mr     <- data_valuation %>% filter(!is.na(SP_Mean_Reversion_Delta)) %>% pull(SP_Mean_Reversion_Delta) %>% sd()
  317.  
  318. ################################################################################
  319. ### Assemble graphs...
  320. ################################################################################
  321.  
  322. g_shiller_pe_z <-
  323.   ggplot(data = data_valuation) +
  324.   theme_bw() +
  325.   geom_line(aes(x = Date, y = CAPE)) +
  326.  
  327.   geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +
  328.  
  329.   geom_hline(yintercept = mean_shiller_pe + 3 * sd_shiller_pe, linetype = "dotted", color = "darkred") +
  330.   geom_hline(yintercept = mean_shiller_pe + 2 * sd_shiller_pe, linetype = "dotted", color = "red") +
  331.   geom_hline(yintercept = mean_shiller_pe + 1 * sd_shiller_pe, linetype = "dotted", color = "firebrick1") +
  332.   geom_hline(yintercept = mean_shiller_pe + 0 * sd_shiller_pe, linetype = "dashed", color = "black") +
  333.   geom_hline(yintercept = mean_shiller_pe - 1 * sd_shiller_pe, linetype = "dotted", color = "blue") +
  334.   geom_hline(yintercept = mean_shiller_pe - 2 * sd_shiller_pe, linetype = "dotted", color = "darkblue") +
  335.  
  336.   labs(x = "",
  337.        y = "",
  338.        #y = "Shiller PE",
  339.        title = "Shiller PE / CAPE") +
  340.   scale_x_date(breaks = pretty_breaks(9)) +
  341.   scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
  342.                      sec.axis = sec_axis(~ (. - mean_shiller_pe) / sd_shiller_pe,
  343.                                          name = "Z-Score",
  344.                                          breaks = pretty_breaks(6)))
  345. print(g_shiller_pe_z)
  346.  
  347.  
  348. g_sp500_mr <-
  349.   ggplot(data = data_valuation) +
  350.   theme_bw() +
  351.   geom_line(aes(x = Date, y = SP_Mean_Reversion_Delta)) +
  352.  
  353.   geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +
  354.  
  355.   geom_hline(yintercept = mean_sp500_mr + 3 * sd_sp500_mr, linetype = "dotted", color = "darkred") +
  356.   geom_hline(yintercept = mean_sp500_mr + 2 * sd_sp500_mr, linetype = "dotted", color = "red") +
  357.   geom_hline(yintercept = mean_sp500_mr + 1 * sd_sp500_mr, linetype = "dotted", color = "firebrick1") +
  358.   geom_hline(yintercept = mean_sp500_mr + 0 * sd_sp500_mr, linetype = "dashed", color = "black") +
  359.   geom_hline(yintercept = mean_sp500_mr - 1 * sd_sp500_mr, linetype = "dotted", color = "blue") +
  360.   geom_hline(yintercept = mean_sp500_mr - 2 * sd_sp500_mr, linetype = "dotted", color = "darkblue") +
  361.  
  362.   labs(x = "",
  363.        y = "",
  364.        title = "Detrended log(Real S&P 500 Price)") +
  365.   scale_x_date(breaks = pretty_breaks(9)) +
  366.   scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
  367.     sec.axis = sec_axis(~ (. - mean_sp500_mr) / sd_sp500_mr,
  368.                         name = "Z-Score",
  369.                         breaks = pretty_breaks(6)))
  370. print(g_sp500_mr)
  371.  
  372.  
  373. g_tobins_q_z <-
  374.   ggplot(data = data_valuation %>% filter(!is.na(Tobins_Q))) +
  375.   theme_bw() +
  376.   geom_line(aes(x = Date, y = Tobins_Q)) +
  377.  
  378.   geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +
  379.  
  380.   geom_hline(yintercept = mean_tobins_q + 3 * sd_tobins_q, linetype = "dotted", color = "darkred") +
  381.   geom_hline(yintercept = mean_tobins_q + 2 * sd_tobins_q, linetype = "dotted", color = "red") +
  382.   geom_hline(yintercept = mean_tobins_q + 1 * sd_tobins_q, linetype = "dotted", color = "firebrick1") +
  383.   geom_hline(yintercept = mean_tobins_q + 0 * sd_tobins_q, linetype = "dashed", color = "black") +
  384.   geom_hline(yintercept = mean_tobins_q - 1 * sd_tobins_q, linetype = "dotted", color = "blue") +
  385.   geom_hline(yintercept = mean_tobins_q - 2 * sd_tobins_q, linetype = "dotted", color = "darkblue") +
  386.  
  387.   labs(x = "",
  388.        y = "",
  389.        #y = "Tobin's Q",
  390.        title = "Tobin's Q") +
  391.   scale_x_date(breaks = pretty_breaks(9)) +
  392.   scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
  393.     sec.axis = sec_axis(~ (. - mean_tobins_q) / sd_tobins_q,
  394.                         name = "Z-Score",
  395.                         breaks = pretty_breaks(6)))
  396. print(g_tobins_q_z)
  397.  
  398.  
  399. g_aiea_z <-
  400.   ggplot(data = data_valuation %>% filter(!is.na(AIEA))) +
  401.   theme_bw() +
  402.   geom_line(aes(x = Date, y = AIEA)) +
  403.  
  404.   geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +
  405.  
  406.   geom_hline(yintercept = mean_aiea + 2 * sd_aiea, linetype = "dotted", color = "red") +
  407.   geom_hline(yintercept = mean_aiea + 1 * sd_aiea, linetype = "dotted", color = "firebrick1") +
  408.   geom_hline(yintercept = mean_aiea + 0 * sd_aiea, linetype = "dashed", color = "black") +
  409.   geom_hline(yintercept = mean_aiea - 1 * sd_aiea, linetype = "dotted", color = "blue") +
  410.   geom_hline(yintercept = mean_aiea - 2 * sd_aiea, linetype = "dotted", color = "darkblue") +
  411.  
  412.   labs(x = "",
  413.        y = "",
  414.        #y = "AIEA",
  415.        title = "Aggregate Investor Equity Allocation") +
  416.   scale_x_date(breaks = pretty_breaks(9)) +
  417.   scale_y_continuous(labels = scales::percent_format(accuracy = 1),
  418.     sec.axis = sec_axis(~ (. - mean_aiea) / sd_aiea,
  419.                         name = "Z-Score",
  420.                         breaks = pretty_breaks(6)))
  421. print(g_aiea_z)
  422.  
  423. g_buffett_z <-
  424.   ggplot(data = data_valuation %>% filter(!is.na(Buffett_Indicator))) +
  425.   theme_bw() +
  426.   geom_line(aes(x = Date, y = Buffett_Indicator)) +
  427.  
  428.   geom_recession_bars(min(data_valuation$Date), max(data_valuation$Date)) +
  429.  
  430.   geom_hline(yintercept = mean_buffett + 3 * sd_buffett, linetype = "dotted", color = "darkred") +
  431.   geom_hline(yintercept = mean_buffett + 2 * sd_buffett, linetype = "dotted", color = "red") +
  432.   geom_hline(yintercept = mean_buffett + 1 * sd_buffett, linetype = "dotted", color = "firebrick1") +
  433.   geom_hline(yintercept = mean_buffett + 0 * sd_buffett, linetype = "dashed", color = "black") +
  434.   geom_hline(yintercept = mean_buffett - 1 * sd_buffett, linetype = "dotted", color = "blue") +
  435.   geom_hline(yintercept = mean_buffett - 2 * sd_buffett, linetype = "dotted", color = "darkblue") +
  436.  
  437.   labs(x = "",
  438.        y = "",
  439.        title = "Buffett Indicator (Total Stock Market Cap / [GDP + Total Fed Assets])") +
  440.   scale_x_date(breaks = pretty_breaks(9)) +
  441.   scale_y_continuous(#labels = scales::percent_format(accuracy = 1),
  442.                      sec.axis = sec_axis(~ (. - mean_buffett) / sd_buffett,
  443.                                          name = "Z-Score",
  444.                                          breaks = pretty_breaks(6)))
  445. print(g_buffett_z)
  446.  
  447. print(
  448.   plot_grid(
  449.     g_shiller_pe_z,
  450.     g_tobins_q_z,
  451.     g_aiea_z,
  452.     g_buffett_z,
  453.     g_sp500_mr,
  454.     nrow = 5, ncol = 1, align = "hv"))
Advertisement
RAW Paste Data Copied
Advertisement