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