################################################################################ ### ### 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"))