Estimating total asset bubble / stock bubble with R

################################################################################
### Derive the size of asset bubbles in the US economy in US 2019 dollars.
###
###    * Mathew Binkley <Mathew.Binkley@Vanderbilt.edu>
###
### Right now this script can derive the total size of the asset bubble,
### and also the size of the stock bubble portion.   I am working on creating
### measures for housing & commercial real estate and commodity bubbles as well.
################################################################################


################################################################################
### Load necessary packages, installing if necessary
################################################################################
packages <- c("fredr", "ggplot2", "ggthemes", "lubridate", "forecast",
              "tsibble", "dplyr")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, quiet = TRUE)
invisible(lapply(packages, "library",
                 quietly = TRUE,
                 character.only = TRUE,
                 warn.conflicts = FALSE))

### Set my FRED API key to access the FRED database.
### You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
fredr_set_key(API_KEY_FRED)

################################################################################
### A data.frame containing recession start/stop dates, so we can add recession
### bars to our graph
################################################################################

recessions.df = read.table(textConnection(
  "Peak, Trough
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"),
  sep=',', colClasses=c('Date', 'Date'), header=TRUE)

### What date ranges do we want?
Date.start = "1996-01-01"
Date.end   = "2017-01-01"  # I wish market cap dataset had more current data...

### FRED Real GDP gives Real GDP in 2012 dollars.  To get it in real
### 2019 dollars, multiply by 1.12 to account for inflation.  Adjust
### inflation if you're running it in different years.
Inflation = 1.12

################################################################################
### Import Household Wealth from FRED
################################################################################
Series = "TNWBSHNO"

data <- as_tsibble(fredr(series_id = Series,
                         frequency = "a",
                         observation_start = as.Date(Date.start),
                         observation_end   = as.Date(Date.end)), index = "date")

Date            <- data %>% pull('date')
HouseholdWealth <- data %>% pull('value')


################################################################################
### Import Stock Market Capitalization to GDP for US from FRED
################################################################################
Series = "DDDM01USA156NWDB"

data <- as_tsibble(fredr(series_id = Series,
                         frequency = "a",
                         observation_start = as.Date(Date.start),
                         observation_end   = as.Date(Date.end)), index = "date")

StockMarketCapitalization_to_GDP <- data %>% pull('value')

################################################################################
### Import Nominal GDP from FRED
################################################################################
Series = "GDP"

data <- as_tsibble(fredr(series_id = Series,
                         frequency = "a",
                         observation_start = as.Date(Date.start),
                         observation_end   = as.Date(Date.end)), index = "date")

NominalGDP <- data %>% pull('value')

################################################################################
### Import Real GDP from FRED
################################################################################
Series = "GDPC1"

data <- as_tsibble(fredr(series_id = Series,
                         frequency = "a",
                         observation_start = as.Date(Date.start),
                         observation_end   = as.Date(Date.end)), index = "date")

RealGDP <- data %>% pull('value')



################################################################################
### Calculation section for total asset bubble size
### (including stocks, homes, gold, etc.)
################################################################################

### Compute the ratio of Household Wealth to Nominal GDP
Wealth_to_GDP    = HouseholdWealth   / NominalGDP

### If you plot Wealth_to_GDP, you see there are some long-term trends that sit
### beneath the asset bubbles.  Judging by their long-term correlation with
### bond prices, they may reflect "bond bubbles", or perhaps longer-term
### economic/public policy changes.
###
### I'm primarily interested in the "frothy" bubbles on top, as they tend not
### to fall below the trends when they pop.   So I extract them from the
### underlying "bond bubbles" by modeling the trend as three piecewise sections,
### taking the local minima inside each trend, and using a linear regression to
### determine the floors.

### Pre-1960:         Wealth_Floor = -72.60170 + 0.038941*Year
### 1960-1978.75:     Wealth_Floor = +52.16769 - 0.024760*Year
### 1982 -> Present:  Wealth_Floor = -63.48539 + 0.033682*Year

AllAssetBubbles <- rep( NA, length( Wealth_to_GDP ))
early <- Date <= as.Date("1959-09-03")
late  <- Date >  as.Date("1978-12-14")
mid   <- !early & !late
AllAssetBubbles[early] <- Wealth_to_GDP[early] + 72.6017 - 0.038941*decimal_date(Date[early])
AllAssetBubbles[mid]   <- Wealth_to_GDP[mid]   - 52.1677 + 0.024760*decimal_date(Date[mid])
AllAssetBubbles[late]  <- Wealth_to_GDP[late]  + 63.4854 - 0.033682*decimal_date(Date[late])

### Calulate the real 2019 dollar value of the asset bubbles
AllAssetBubbles_RealDollars = Inflation * AllAssetBubbles * RealGDP / 1000

################################################################################
### Calculation section for stock bubble size
################################################################################

### Analogous to extracting total asset bubble size by examining household
### wealth to GDP, we examine stock market wealth to GDP to find the size of the
### stock market bubble.  This only approximate for several reasons such as
### significantly less data to analyze with, the coarseness of the data
### (annually for household wealth vs quarterly for stock market cap from FRED),
### and the initial starting date of the value is close to the dot.com bubble.
### Better data sources would help tremendously here.

### Get the Stock Market Cap to GDP, scaled correctly for our calculation
StockBubble_in_GDP = ((StockMarketCapitalization_to_GDP/100) - 1)

### Multiply by Real GDP, multiply by inflation, and scale by 1000
### to get the size of the stock bubble in trillions of real 2019 US $
StockBubbleSize = Inflation * StockBubble_in_GDP * RealGDP / 1000


################################################################################
### Calculation section for housing bubble size
################################################################################

### Note:  I'm still working on a method to directly derive the size of the
### housing bubble.   Until then, as a decent approximation, assume that:
HousingBubble = AllBubbles - StockBubbles


################################################################################
### Graphing section
################################################################################

### Our data frame for graphing...
data.df = data.frame(Date = Date,
                     AllBubbles = AllAssetBubbles_RealDollars,
                     StockBubbles = StockBubbleSize,
                     HousingBubbles = HousingBubble)

### The subset of recessions that lie within the date range of our data
recessions.trim = subset(recessions.df, Peak >= min(Date))

### Graph globals
AnnotationColor = "black"

### Common caption strings to make doing the caption easier
C1 = "Board of Governors of the Federal Reserve System (US), Households and nonprofit organizations; net worth, Level [TNWBSHNO]\n"
C2 = "World Bank, Stock Market Capitalization to GDP for United States [DDDM01USA156NWDB]\n"
C3 = "U.S. Bureau of Economic Analysis, Gross Domestic Product [GDP]\n"
C4 = "U.S. Bureau of Economic Analysis, Real Gross Domestic Product [GDPC1]\n"
C5 = paste("Data retrieved from FRED, Federal Research Bank of St. Louis on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))

### Create the graph
Title    = "Total Asset Bubble Size alongside Stock Bubble Size"
Subtitle = "Recessions marked with vertical bars"
Caption  = paste(C1, C2, C3, C4, C5)
XLab     = "Year"
YLab     = "Trillions of 2019 US $"

P <- ggplot(data = data.df, mapping = aes(x = Date, y = AllBubbles)) +

   theme_economist() + scale_colour_economist() +
   theme(legend.title = element_blank()) +
   labs(title = Title, subtitle = Subtitle, caption = Caption, x = XLab, y = YLab) +
   scale_x_date(limits = c(as.Date(Date.start), as.Date(Date.end))) +


   geom_line(data = data.df, size = 1.3,
             aes(y = AllBubbles,
                 color = "All Asset Bubbles", linetype = "All Asset Bubbles")) +

   geom_line(data = data.df, size = 1.3,
             aes(y = StockBubbles,
                 color = "Stock Bubbles", linetype = "Stock Bubbles")) +

   geom_line(data = data.df, size = 1.3,
             aes(y = AllBubbles - StockBubbles,
                 color = "Imputed Housing Bubbles", linetype = "Imputed Housing Bubbles")) +

   scale_linetype_manual(name = "colour",
                         values = c("All Asset Bubbles"       = "solid",
                                    "Stock Bubbles"           = "solid",
                                    "Imputed Housing Bubbles" = "dotted")) +

   geom_rect(data = recessions.trim, inherit.aes = F, fill='darkgray', alpha=0.25,
             aes(xmin = as.Date(Peak), xmax = as.Date(Trough), ymin = -Inf,  ymax = +Inf)) +

   annotate("text", x = as.Date("2000-01-01"), y = 10,
            label = "Dot.Com\nBubble", size  = 5, color = AnnotationColor) +
   annotate("text", x = as.Date("2006-04-01"), y = 16,
            label = "Great Recession\nBubble", size  = 5, color = AnnotationColor) +
   annotate("text", x = as.Date("2017-01-01"), y = 19,
            label = "\nCurrent\nBubble", size  = 5, color = AnnotationColor)
print(P)