Stocks, Bonds, Bills, and Inflation’s Components: Online Appendix

Libraries and theme

Loads libraries and custom theme

library(tidyquant)

## Loading required package: lubridate

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## Loading required package: PerformanceAnalytics

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

## Loading required package: quantmod

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## == Need to Learn tidyquant? ====================================================
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x dplyr::filter()          masks stats::filter()
## x dplyr::first()           masks xts::first()
## x lubridate::intersect()   masks base::intersect()
## x dplyr::lag()             masks stats::lag()
## x dplyr::last()            masks xts::last()
## x lubridate::setdiff()     masks base::setdiff()
## x lubridate::union()       masks base::union()

library(broom)
library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(tsibble)

## 
## Attaching package: 'tsibble'

## The following object is masked from 'package:zoo':
## 
##     index

## The following object is masked from 'package:lubridate':
## 
##     interval

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

library(ggcorrplot)
library(ggfittext)
library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(data.table)

## 
## Attaching package: 'data.table'

## The following object is masked from 'package:tsibble':
## 
##     key

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following objects are masked from 'package:xts':
## 
##     first, last

## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year

library(roll)
library(ivreg)
library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(readxl)
library(fable)

## Loading required package: fabletools

## 
## Attaching package: 'fable'

## The following object is masked from 'package:tidyquant':
## 
##     VAR

library(mFilter)
library(ggrepel)
library(lmtest)
library(magrittr)

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

library(sandwich)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

library(slider)
library(tsibble)
library(feasts)
library(tseries)
library(ggfittext) # for geom_fit_text
library(ggalt) # for geom_encircle

## Registered S3 methods overwritten by 'ggalt':
##   method                  from   
##   grid.draw.absoluteGrob  ggplot2
##   grobHeight.absoluteGrob ggplot2
##   grobWidth.absoluteGrob  ggplot2
##   grobX.absoluteGrob      ggplot2
##   grobY.absoluteGrob      ggplot2

library(stargazer) # for markdown printing: {r results = "asis", message=FALSE, warning=FALSE}

## 
## Please cite as:

##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

library(fixest)
library(patchwork)


theme =  
  theme(
    plot.title = element_text(size = 30, lineheight = .9, face = "bold", hjust = .5, vjust = .75, family = "serif"
    ),
    plot.subtitle = element_text(hjust = .5, vjust = 1, family = "serif", size = 20),
    plot.caption = element_text(hjust = .5, family = "serif", size = 15),
    panel.background = element_rect(fill = "white"),
    axis.title = element_text(family = "serif", size = 20),
    axis.text = element_text(size = 20, family = "serif"),
    axis.line.y = element_line(color = "black"),
    axis.line.x = element_line(color = "black"),
    strip.text = element_text(size = 15, family = "serif", face = "bold"),
    strip.background.x = element_rect(fill = "lightblue"))

One-year data

Supply chain disruptions

Monthly data

sbbi_mo <- 
  read_csv("sbbi.csv", locale = locale(encoding = "latin1")) %>% 
    select(c(1, 2, 5, 7, 15)) %>% 
    set_colnames(c("date", "large", "small", "lt_gov", "tbills")) %>% 
  mutate(date = seq(as.Date("1926-01-01"),
                    by = "month",
                    length = length(large))) 

cpi_short_mo <- 
  c("CPIAUCSL", "CPILFESL") %>% 
  tq_get(get = "economic.data", from = "1957-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  select(date, cpi = cpiaucsl, cpi_core = cpilfesl) %>% 
  mutate(across(2:3, ~ (./lag(., 1)) - 1)) %>% 
  drop_na 

median_cpi_mo <- 
  c("MEDCPIM158SFRBCLE") %>% 
  tq_get(get = "economic.data", from = "1983-01-01") %>% 
  select(date, median = price) %>%
  mutate(median_cpi = ((1 + (median / 100))^(1/12)) - 1) %>% 
  drop_na %>% 
  select(date, median_cpi)

trimmed_cpi_mo <- 
  c("TRMMEANCPIM158SFRBCLE") %>% 
  tq_get(get = "economic.data", from = "1983-01-01") %>% 
  select(date, trim_cpi = price) %>% 
  mutate(trimmed_cpi = ((1 + (trim_cpi / 100)) ^ (1/12)) - 1) %>% 
  select(date, trimmed_cpi) %>% 
  drop_na 

pce_mo <- 
  c("PCEPI", "PCEPILFE") %>% 
  tq_get(get = "economic.data", from = "1957-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  mutate(across(2:3, ~ (./lag(., 1)) - 1)) %>% 
  drop_na 

trimmed_pce_mo <- 
  c("PCETRIM12M159SFRBDAL") %>% 
  tq_get(get = "economic.data", from = "1978-01-01") %>% 
  select(date, trimmed_pce = price) %>% 
  mutate(trimmed_pce =  ((1 + trimmed_pce/100)^(1/12)) - 1)

median_pce_mo <- 
  read_csv("median_pce.csv") %>% 
  mutate(date = mdy(date)) %>% 
  drop_na %>% 
  select(date, median_pce = med_pce) %>% 
  mutate(median_pce = ((1 + median_pce/100)^(1/12)) - 1)


russell1 <- 
  read_csv("russell.csv") %>% 
    mutate(russell1 = russell1/100) %>% 
  mutate(date = floor_date(mdy(date), "month")) %>% 
  arrange(date)

russell2 <- 
  read_csv("russell2.csv") %>% 
    mutate(russell2 = russell2/100) %>% 
  mutate(date = floor_date(mdy(date), "month")) %>% 
  arrange(date)

  tightening <- 
    c("FEDFUNDS") %>% 
    tq_get(get = "economic.data", from = "1954-01-01") %>% 
    select(date, fedfunds = price) %>%
    mutate(tighteninga = ifelse(fedfunds >= lag(fedfunds, 6) + .5, 1, 0),
           tighteningb = ifelse(lag(tighteninga ,1) == 1 & fedfunds >= lag(fedfunds) - .15, 1, 0),
           tightening = factor(ifelse(tighteninga == 1 | tighteningb == 1, 1, 0)),
           year = year(date)) %>% 
    group_by(year, tightening) %>% 
    mutate(first = first(date), last = last(date)) %>% 
    ungroup %>% 
    select(date, tightening) %>% 
    drop_na
    #filter(tightening == 1)
  
  # BE infl
  
  real <- 
   c("REAINTRATREARAT10Y") %>% 
   tq_get(get = "economic.data", from = "1947-01-01") %>% 
   select(date, real_yield = price)
 
 ten <- 
   c("DGS10") %>% 
   tq_get(get = "economic.data", from = "1947-01-01") %>% 
   select(date, yield = price) %>%
   mutate(date = floor_date(date, "month")) %>% 
   fill(yield, .direction = c("up")) %>% 
   group_by(date) %>% 
   filter(row_number() == n()) %>% 
   ungroup
   
   be_infl <- 
   real %>% 
   left_join(ten, by = "date") %>% 
   mutate(be_infl = yield - real_yield) %>% 
   select(date, be_infl)
   
     tcu <- 
    c("TCU") %>% 
    tq_get(get = "economic.data", from = "1947-01-01") %>% 
    select(date, tcu = price) %>% 
    mutate(hot_tcu = factor(ifelse(tcu > 82, 1, 0))) %>% 
    drop_na %>% 
    select(date, hot_tcu)
     
      unrate <- 
   c("UNRATE", "NROU") %>% 
   tq_get(get = "economic.data", from = "1949-01-01") %>% 
   pivot_wider(names_from = symbol, values_from = price) %>% 
   clean_names %>% 
   fill(nrou, .direction = "up") %>% 
   mutate(hot_labor = factor(ifelse(unrate <= nrou, 1, 0))) %>% 
   select(date, hot_labor)
      
    ismmsdi <- 
    read_csv("ismmsdi.csv") %>% 
      mutate(date = floor_date(mdy(date), "month")) %>% 
   select(date, ismmsdi) %>% 
      mutate(ismmsdi_diff = ismmsdi - lag(ismmsdi)) %>% 
      select(date, ismmsdi_diff)
    
    
supply <- 
  read_csv("supply.csv") %>% 
  mutate(date = floor_date(dmy(date), "month"),
         gscpi_diff = gscpi - lag(gscpi)) %>% 
  select(date, gscpi_diff)

 
 mich <- 
   c("MICH") %>% 
  tq_get(get = "economic.data", from = "1977-01-01") %>% 
  select(date, mich = price) 
 
 cleve <- 
 c("EXPINF1YR") %>% 
    tq_get(get = "economic.data", from = "1977-01-01") %>% 
  select(date, cleve = price) 
 
 # real PCE
 rpcets <- 
  c("PCE", "CPIAUCSL") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  mutate(rpce = 100*pce/cpiaucsl) %>% 
  select(rpce) %>% 
  ts(start = 1959, frequency = 12)
 
  filter <- 
  na.omit(rpcets) %>% 
  hpfilter(drift = F, type = c("lambda"), freq = 129600) 
  
  rpce <- 
    tibble(trend = as.vector(filter$trend),
       cycle = as.vector(filter$cycle),
       date = seq(as.Date("1959-01-01"),
                              by = "month", 
                              length = length(filter$trend))) %>% 
    mutate(month = floor_date(date, "month")) %>% 
    select(date = month, rpce_cycle = cycle) %>% 
    select(date, rpce_cycle)
  
  # employment level
  
  payrolls <- 
  c("PAYEMS") %>% 
  tq_get(get = "economic.data", from = "1939-01-01") %>% 
  select(empl = price) %>% 
  ts(start = 1939, frequency = 12)
 
  filter <- 
  payrolls %>% 
  hpfilter(drift = F, type = c("lambda"), freq = 129600) 
  
  payrolls <- 
    tibble(trend = as.vector(filter$trend),
       cycle = as.vector(filter$cycle),
       date = seq(as.Date("1939-01-01"),
                              by = "month", 
                              length = length(filter$trend))) %>% 
    mutate(month = floor_date(date, "month")) %>% 
    select(date = month, payrolls_cycle = cycle) %>% 
    select(date, payrolls_cycle)
  
  
  # CPI stochastic trend and cycle
  
  pce <- 
    c("PCEPI") %>% 
    tq_get(get = "economic.data", from = "1949-01-01") %>% 
    select(date, pce = price) %>% 
    mutate(date = yearmonth(date)) %>% 
    as_tsibble
  
  trend <- 
    pce %>% 
    model(stochastic = ARIMA(log(pce) ~ pdq(d = 1)))
  
  pce_def_cyc <- 
    trend %>% 
    augment %>% 
    select(date, pce, .fitted, .resid, .innov) %>%
    mutate(cycle = .resid + .innov) %>% 
    select(date, pce_def_cycle = cycle, pce_def_trend = .fitted) %>% 
    mutate(date = ym(date),
           pce_def_trend = pce_def_trend - lag(pce_def_trend, 1))
  

df_common_mo <- 
sbbi_mo %>% 
  right_join(russell1, by = "date") %>% 
  right_join(russell2, by = "date") %>% 
  left_join(cpi_short_mo, by = "date") %>% 
  right_join(median_cpi_mo, by = "date") %>% 
  right_join(trimmed_cpi_mo, by = "date") %>% 
  right_join(pce_mo, by = "date") %>% 
  right_join(median_pce_mo, by = "date") %>% 
  right_join(trimmed_pce_mo, by = "date") %>% 
    mutate(total_median_shock = cpi - median_cpi, # all rel price changes
           total_trimmed_mean_shock = cpi - trimmed_cpi,# all rel price changes
           fe_shock = cpi - cpi_core,
           other_median_shock = cpi_core - median_cpi) %>% 
  left_join(tightening, by = "date") %>% 
  left_join(be_infl, by = "date") %>% 
  left_join(tcu, by = "date") %>% 
  left_join(unrate, by = "date") %>% 
  left_join(mich, by = "date") %>% 
  left_join(cleve, by = "date") %>% 
  left_join(ismmsdi, by = "date") %>% 
  left_join(supply, by = "date") %>% 
  left_join(rpce, by = "date") %>% 
  left_join(pce_def_cyc, by = "date") %>% 
  mutate(time = row_number()) %>% 
  filter(year(date) < 2023)

Figure 1: large cap stocks and inflation

Figure 2 stocks and headline vs. median inflation

plot1 <- 
  df_common_mo %>% 
  drop_na %>% 
  filter(year(date) > 1999) %>% 
  ggplot(aes(x = median_cpi, y = russell1)) +
  geom_point() +
  geom_smooth(method = "lm", color = "darkgray", se = F) +
  theme +
  labs(x = "Median CPI", y = "Large cap return")

plot2 <- 
  df_common_mo %>% 
  filter(year(date) > 1999) %>% 
  drop_na %>% 
  ggplot(aes(x = total_median_shock, y = russell1)) +
  geom_point() +
  geom_smooth(method = "lm", color = "darkgray", se = F) +
  theme +
  labs(x = "Headline CPI shock", y = "Large cap return")

plot1 + plot2 + plot_annotation(caption = "Source: CFA Institute, FRED")

## Table 1 Summary statistics (for first order autocorrelations)

  data <- 
    df_common_mo %>% #make sure categorical vars are excluded
    select(2:19, rpce_cycle, pce_def_cycle, pce_def_trend) %>% 
  summarize(count = across(where(is.numeric), ~ n()),
            mean = across(where(is.numeric), ~ mean(., na.rm = T)),
            sd = across(where(is.numeric), ~ sd(., na.rm = T)),
            iqr = across(where(is.numeric), ~ IQR(., na.rm = T)),
            autocor = across(where(is.numeric), ~ cor(x = ., y = lag(., 1), 
                                                      use = "complete.obs"))) %>% 
  mutate(index = 1) %>% 
  relocate(index) %>% 
  pivot_longer(-index) %>% 
  select(-index) 
             
  tibble(statistic = c("obs", "mean", "sd", "IQR", "autocorr"), data$value[,1:18]) %>% 
    pivot_longer(cols = large:other_median_shock) %>% 
    pivot_wider(names_from = statistic, values_from = value) %>% 
    rename(series = name) %>% 
    kbl(digits = 3) %>% 
    kable_classic(full_width = F)

series	obs	mean	sd	IQR	autocorr
large	540	0.010	0.044	0.052	0.003
small	540	0.010	0.058	0.066	0.117
lt_gov	540	0.007	0.031	0.039	0.068
tbills	540	0.003	0.002	0.004	0.983
russell1	540	0.010	0.044	0.054	0.019
russell2	540	0.009	0.056	0.068	0.068
cpi	540	0.002	0.003	0.002	0.478
cpi_core	540	0.002	0.001	0.002	0.581
median_cpi	540	0.002	0.001	0.001	0.726
trimmed_cpi	540	0.002	0.001	0.001	0.706
pcepi	540	0.002	0.002	0.002	0.665
pcepilfe	540	0.002	0.002	0.002	0.710
median_pce	540	0.003	0.001	0.001	0.998
trimmed_pce	540	0.002	0.001	0.001	0.998
total_median_shock	540	0.000	0.002	0.002	0.423
total_trimmed_mean_shock	540	0.000	0.002	0.001	0.400
fe_shock	540	0.000	0.002	0.002	0.360
other_median_shock	540	0.000	0.001	0.001	0.288

Figure appendix: Fed tightenings

tightening_shading <- 
  c("FEDFUNDS") %>% 
    tq_get(get = "economic.data", from = "1954-01-01") %>% 
    select(date, fedfunds = price) %>%
    mutate(tighteninga = ifelse(fedfunds >= lag(fedfunds, 6) + .5, 1, 0),
           tighteningb = ifelse(lag(tighteninga ,1) == 1 & fedfunds >= lag(fedfunds) - .15, 1, 0),
           tightening = ifelse(tighteninga == 1 | tighteningb == 1, 1, 0),
           year = year(date)) %>% 
    group_by(year, tightening) %>% 
    mutate(first = first(date), last = last(date)) %>% 
    ungroup %>% 
    filter(year(date) > 1982) %>% 
    filter(tightening == 1) 


fedfunds %>% 
    ggplot(aes(x = date, y = fedfunds)) +
    geom_rect(aes(xmin = first, xmax = last, ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_rect(aes(xmin = as.Date("1988-06-01"), xmax = as.Date("1989-06-01"), ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_rect(aes(xmin = as.Date("1994-04-01"), xmax = as.Date("1995-07-01"), ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_rect(aes(xmin = as.Date("1999-11-01"), xmax = as.Date("2000-10-01"), ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_rect(aes(xmin = as.Date("2004-09-01"), xmax = as.Date("2006-10-01"), ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_rect(aes(xmin = as.Date("2022-05-01"), xmax = as.Date("2023-11-01"), ymin = -Inf, ymax = +Inf), 
              data = tightening_shading, fill = "pink") +
    geom_line() +
    theme +
    labs(x = NULL, y = "Fed funds rate", caption = "Source: FRED, Author calculations")

Figure 3: Stocks and inflation during fed tightenings

df_common_mo %>% 
  ggplot(aes(x = cpi, y = large)) +
  geom_point(aes(shape = factor(tightening)), size = 3) +
  geom_smooth(aes(linetype = factor(tightening)), method = "lm", se = F, color = "black", show.legend = T) +
  scale_shape_manual(values = c(1, 3), labels = c("N", "Y")) +
  scale_linetype_manual(values = c(1, 3), labels = c("N", "Y")) +
  theme +
  labs(x = "CPI one-month change", y = "Large cap one-month change", 
       caption = "Source: CFA Institute, FRED") +
  guides(shape = guide_legend(title = "Fed tightening?"),
         linetype = guide_legend(title = "Fed tightening?"))

Figure NA: Correlation plot

# 1983-2022
 
 df_common_mo %>% 
  select(-date, -tightening, -hot_tcu, -hot_labor) %>%
  cor(use = "complete.obs") %>% 
   ggcorrplot(type = "upper", show.legend = T, lab = T, colors = c("gray", "lightgray", "white")) +
   theme 
  #round(2)
  #stargazer(type = "text", out = "corr.html", digits = 2)

Regressions: OLS using Ibbotson

# CPI

mod_large_cpi <- 
  df_common_mo %>% 
lm(large ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
lm(small ~ cpi, data = .) 

mod_ltbond_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
lm(tbills ~ cpi, data = .) 

# Core CPI

mod_large_cpi_core <- 
  df_common_mo %>% 
lm(large ~ cpi_core, data = .) 

mod_small_cpi_core <- 
 df_common_mo %>% 
lm(small ~ cpi_core, data = .) 

mod_ltbond_cpi_core <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi_core, data = .)

mod_tbills_cpi_core <- 
  df_common_mo %>% 
lm(tbills ~ cpi_core, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
lm(large ~ median_cpi, data = .) 

mod_small_median_cpi<- 
 df_common_mo %>% 
lm(small ~ median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
lm(large ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
lm(small ~ median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
lm(large ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
lm(small ~ trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  lm(large ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  lm(small ~ total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
lm(large ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
lm(small ~ pce_def_cycle, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_cycle, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
lm(large ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
lm(small ~ pce_def_trend, data = .) 

mod_ltbond_pce_def_trend <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_trend, data = .) 

mod_tbills_pce_def_trend <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
       sqrt(diag(vcovHC(mod_large_cpi_core))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_cpi_core))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi))),
              sqrt(diag(vcovHC(mod_ltbond_cpi_core))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
               sqrt(diag(vcovHAC(mod_tbills_cpi_core))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_cpi_core,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_cpi_core,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_cpi_core,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_cpi_core,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "HC", "IID", "HC", "IID", "HAC", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_OLS_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap								Small cap								LT gov. bond				T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)	(21)	(22)	(23)	(24)
cpi	-0.149								0.264								-2.759***				0.153
	(1.171)								(1.584)								(0.777)				(0.094)
cpi_core		-1.271								-3.449								-0.475				0.729**
		(1.801)								(2.181)								(1.042)				(0.340)
median_cpi			-3.265*								-5.699**								-1.731				1.054*
			(1.936)								(2.536)								(1.356)				(0.549)
median_pce				-2.279								-4.609
				(3.113)								(3.547)
trimmed_cpi					-1.974								-3.760*
					(1.731)								(2.270)
total_median_shock						0.418								1.367						-3.036***				-0.008
						(1.263)								(1.701)						(0.565)				(0.062)
pce_def_trend							-0.010								-0.015
							(0.009)								(0.014)
pce_def_cycle								0.023								0.045
								(0.022)								(0.028)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	HC	IID	IID	IID	HC	HC	HC	HC	HC	IID	HC	IID	HAC	HAC	HAC	HAC
Observations	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480
R2	0.0001	0.002	0.006	0.001	0.003	0.001	0.002	0.006	0.0001	0.007	0.010	0.004	0.006	0.003	0.003	0.014	0.057	0.0005	0.003	0.057	0.030	0.194	0.216	0.0001
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 1983-1999


mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ cpi, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ cpi, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_cpi, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(small ~ median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(large ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(small ~ total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_cycle, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ pce_def_cycle, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(tbills ~ pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_1983_1999_OLS_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-3.723**							-4.469**							-2.035			0.235***
	(1.836)							(2.220)							(1.451)			(0.084)
median_cpi		-2.833							-5.727							0.939			1.036***
		(3.467)							(4.182)							(2.334)			(0.153)
median_pce			-2.269							-2.831
			(4.593)							(5.498)
trimmed_cpi				-2.702							-4.077
				(3.057)							(3.695)
total_median_shock					-4.031*							-4.137					-2.953**			-0.004
					(2.057)							(2.519)					(1.266)			(0.071)
pce_def_trend						-0.053**							-0.046
						(0.024)							(0.030)
pce_def_cycle							-0.032							-0.042
							(0.041)							(0.051)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204
R2	0.025	0.003	0.001	0.004	0.022	0.023	0.004	0.024	0.009	0.001	0.006	0.016	0.012	0.005	0.016	0.001	0.026	0.072	0.319	0.00001
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 2000 - 2022



mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ cpi, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ cpi, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_cpi, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(large ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(small ~ total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(small ~ pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_2000_2022_OLS_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	0.467							1.326							-3.083***			0.040
	(1.312)							(1.794)							(0.908)			(0.037)
median_cpi		-6.162**							-7.269**							-4.249**			0.351
		(2.504)							(3.484)							(1.805)			(0.287)
median_pce			-11.050*							-11.202*
			(6.161)							(5.919)
trimmed_cpi				-3.785*							-4.747
				(2.249)							(3.122)
total_median_shock					1.391							2.553					-3.045***			-0.001
					(1.372)							(1.885)					(0.648)			(0.034)
pce_def_trend						-0.002							-0.010
						(0.010)							(0.016)
pce_def_cycle							0.032							0.059*
							(0.024)							(0.031)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276
R2	0.001	0.022	0.024	0.010	0.008	0.0001	0.017	0.004	0.016	0.013	0.008	0.014	0.002	0.030	0.090	0.020	0.075	0.007	0.064	0.00000
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: OLS using Ibbotson with tightening

# CPI

mod_large_cpi <- 
  df_common_mo %>% 
lm(large ~ cpi + tightening, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
lm(small ~ cpi + tightening, data = .) 

mod_ltbond_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi + tightening, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
lm(tbills ~ cpi + tightening, data = .) 

# Core CPI

mod_large_cpi_core <- 
  df_common_mo %>% 
lm(large ~ cpi_core + tightening, data = .) 

mod_small_cpi_core <- 
 df_common_mo %>% 
lm(small ~ cpi_core + tightening, data = .) 

mod_ltbond_cpi_core <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi_core + tightening, data = .)

mod_tbills_cpi_core <- 
  df_common_mo %>% 
lm(tbills ~ cpi_core + tightening, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
lm(large ~ median_cpi + tightening, data = .) 

mod_small_median_cpi<- 
 df_common_mo %>% 
lm(small ~ median_cpi + tightening, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ median_cpi + tightening, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
lm(tbills ~ median_cpi + tightening, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
lm(large ~ median_pce + tightening, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
lm(small ~ median_pce + tightening, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
lm(lt_gov ~ median_pce + tightening, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
lm(large ~ trimmed_cpi + tightening, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
lm(small ~ trimmed_cpi + tightening, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
lm(lt_gov ~ trimmed_cpi + tightening, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  lm(large ~ total_median_shock + tightening, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  lm(small ~ total_median_shock + tightening, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ total_median_shock + tightening, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ total_median_shock + tightening, data = .) 

# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
lm(large ~ pce_def_cycle + tightening, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
lm(small ~ pce_def_cycle + tightening, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_cycle + tightening, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_cycle + tightening, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
lm(large ~ pce_def_trend + tightening, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
lm(small ~ pce_def_trend + tightening, data = .) 

mod_ltbond_pce_def_trend <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_trend + tightening, data = .) 

mod_tbills_pce_def_trend <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_trend + tightening, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
       sqrt(diag(vcovHC(mod_large_cpi_core))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_cpi_core))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi))),
              sqrt(diag(vcovHC(mod_ltbond_cpi_core))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
               sqrt(diag(vcovHAC(mod_tbills_cpi_core))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_cpi_core,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_cpi_core,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_cpi_core,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_cpi_core,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "HC", "IID", "HC", "IID", "HAC", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_OLS_Ibbotson_tightening.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap								Small cap								LT gov. bond				T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)	(21)	(22)	(23)	(24)
cpi	-0.099								0.316								-2.792***				0.127
	(1.172)								(1.585)								(0.778)				(0.092)
cpi_core		-1.151								-3.405								-0.490				0.666**
		(1.771)								(2.188)								(1.026)				(0.293)
median_cpi			-3.157								-5.783**								-1.855				0.948*
			(2.004)								(2.626)								(1.403)				(0.511)
median_pce				-1.977								-4.487
				(3.068)								(3.662)
trimmed_cpi					-1.824								-3.708
					(1.772)								(2.324)
total_median_shock						0.423								1.372						-3.036***				-0.010
						(1.268)								(1.706)						(0.566)				(0.062)
pce_def_trend							-0.009								-0.014
							(0.009)								(0.014)
pce_def_cycle								0.023								0.045
								(0.022)								(0.029)
tightening1	-0.003	-0.002	-0.001	-0.002	-0.002	-0.003	-0.003	-0.003	-0.003	-0.001	0.001	-0.001	-0.001	-0.003	-0.002	-0.002	0.002	0.0003	0.001	0.0002	0.002***	0.001**	0.001**	0.002***
	(0.005)	(0.004)	(0.005)	(0.005)	(0.005)	(0.005)	(0.005)	(0.005)	(0.007)	(0.007)	(0.006)	(0.006)	(0.006)	(0.007)	(0.007)	(0.007)	(0.003)	(0.003)	(0.003)	(0.003)	(0.001)	(0.0005)	(0.0005)	(0.001)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	HC	IID	IID	IID	HC	HC	HC	HC	HC	IID	HC	IID	HAC	HAC	HAC	HAC
Observations	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480
Adjusted R2	-0.003	-0.002	0.002	-0.002	-0.001	-0.003	-0.001	0.003	-0.004	0.003	0.006	-0.001	0.002	-0.0005	-0.001	0.010	0.054	-0.004	-0.001	0.053	0.100	0.238	0.244	0.080
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 1983-1999


mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ cpi + tightening, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ cpi + tightening, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ cpi + tightening, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ cpi + tightening, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_cpi + tightening, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ median_cpi + tightening, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_cpi + tightening, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ median_cpi + tightening, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_pce + tightening, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(small ~ median_pce + tightening, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_pce + tightening, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ trimmed_cpi + tightening, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ trimmed_cpi + tightening, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ trimmed_cpi + tightening, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(large ~ total_median_shock + tightening, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(small ~ total_median_shock + tightening, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ total_median_shock + tightening, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ total_median_shock + tightening, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ pce_def_cycle + tightening, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_cycle + tightening, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ pce_def_cycle + tightening, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(tbills ~ pce_def_cycle + tightening, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ pce_def_trend + tightening, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_trend + tightening, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_1983_1999_OLS_Ibbotson_tightening.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-3.621*							-4.395*							-2.065			0.202**
	(1.914)							(2.329)							(1.466)			(0.091)
median_cpi		-2.393							-5.449							1.065			0.976***
		(3.560)							(4.297)							(2.398)			(0.124)
median_pce			-1.758							-2.352
			(4.509)							(5.589)
trimmed_cpi				-2.361							-3.800
				(3.121)							(3.773)
total_median_shock					-3.952*							-4.061					-2.955**			-0.023
					(2.098)							(2.567)					(1.272)			(0.073)
pce_def_trend						-0.052**							-0.045
						(0.025)							(0.031)
pce_def_cycle							-0.031							-0.041
							(0.041)							(0.052)
tightening1	-0.002	-0.004	-0.005	-0.004	-0.004	-0.003	-0.005	-0.002	-0.002	-0.004	-0.003	-0.004	-0.004	-0.005	0.001	-0.001	0.0001	0.001	0.001	0.001
	(0.007)	(0.007)	(0.007)	(0.007)	(0.007)	(0.007)	(0.007)	(0.009)	(0.008)	(0.008)	(0.008)	(0.009)	(0.009)	(0.009)	(0.005)	(0.005)	(0.005)	(0.001)	(0.0004)	(0.001)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204
Adjusted R2	0.016	-0.005	-0.007	-0.004	0.014	0.015	-0.003	0.015	-0.0002	-0.007	-0.003	0.007	0.003	-0.004	0.007	-0.009	0.017	0.113	0.333	0.062
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 2000 - 2022



mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ cpi + tightening, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ cpi + tightening, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ cpi + tightening, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ cpi + tightening, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_cpi + tightening, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ median_cpi + tightening, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_cpi + tightening, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(tbills ~ median_cpi + tightening, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_pce + tightening, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ median_pce + tightening, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_pce + tightening, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ trimmed_cpi + tightening, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ trimmed_cpi + tightening, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ trimmed_cpi + tightening, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(large ~ total_median_shock + tightening, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(small ~ total_median_shock + tightening, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ total_median_shock + tightening, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ total_median_shock + tightening, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_cycle + tightening, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ pce_def_cycle + tightening, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_trend + tightening, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(small ~ pce_def_trend + tightening, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_2000_2022_OLS_Ibbotson_tightening.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	0.495							1.350							-3.104***			0.022
	(1.309)							(1.786)							(0.909)			(0.042)
median_cpi		-6.397**							-7.705**							-4.619**			0.177
		(2.610)							(3.630)							(1.880)			(0.270)
median_pce			-11.788*							-12.185*
			(6.190)							(6.267)
trimmed_cpi				-3.782							-4.862
				(2.306)							(3.202)
total_median_shock					1.386							2.550					-3.045***			0.003
					(1.381)							(1.897)					(0.649)			(0.032)
pce_def_trend						-0.002							-0.010
						(0.010)							(0.016)
pce_def_cycle							0.032							0.059*
							(0.025)							(0.032)
tightening1	-0.003	0.002	0.004	-0.00005	-0.002	-0.002	-0.002	-0.002	0.004	0.005	0.002	-0.001	-0.001	-0.0002	0.002	0.004	-0.0002	0.002***	0.002***	0.002***
	(0.006)	(0.007)	(0.006)	(0.007)	(0.006)	(0.007)	(0.007)	(0.010)	(0.010)	(0.010)	(0.009)	(0.010)	(0.010)	(0.010)	(0.004)	(0.005)	(0.005)	(0.0005)	(0.0005)	(0.0005)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276
Adjusted R2	-0.006	0.015	0.018	0.003	0.001	-0.007	0.010	-0.003	0.009	0.007	0.001	0.007	-0.005	0.023	0.084	0.014	0.068	0.236	0.249	0.234
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: OLS using Ibbotson with time

# CPI

mod_large_cpi <- 
  df_common_mo %>% 
lm(large ~ cpi + time, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
lm(small ~ cpi + time, data = .) 

mod_ltbond_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi + time, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
lm(tbills ~ cpi + time, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
lm(large ~ median_cpi + time, data = .) 

mod_small_median_cpi<- 
 df_common_mo %>% 
lm(small ~ median_cpi + time, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ median_cpi + time, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
lm(tbills ~ median_cpi + time, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
lm(large ~ median_pce + time, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
lm(small ~median_pce + time, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
lm(lt_gov ~ median_pce + time, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
lm(large ~ trimmed_cpi + time, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
lm(small ~trimmed_cpi + time, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
lm(lt_gov ~ trimmed_cpi + time, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  lm(large ~ total_median_shock + time, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  lm(small ~total_median_shock + time, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ total_median_shock + time, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ total_median_shock + time, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
lm(large ~ pce_def_cycle + time, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
lm(small ~pce_def_cycle + time, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_cycle + time, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_cycle + time, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
lm(large ~ pce_def_trend + time, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
lm(small ~ pce_def_trend + time, data = .) 

mod_ltbond_pce_def_trend <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_trend + time, data = .) 

mod_tbills_pce_def_trend <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_trend + time, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_time_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-0.226							0.224							-2.910***			0.080***
	(1.170)							(1.578)							(0.778)			(0.011)
median_cpi		-4.067**							-6.442**							-2.770**			0.589***
		(2.005)							(2.631)							(1.396)			(0.120)
median_pce			-5.529							-7.954*
			(4.409)							(4.277)
trimmed_cpi				-2.581							-4.284*
				(1.789)							(2.348)
total_median_shock					0.421							1.369					-3.030***			-0.004
					(1.261)							(1.701)					(0.562)			(0.013)
pce_def_trend						-0.009							-0.014
						(0.009)							(0.014)
pce_def_cycle							0.023							0.046
							(0.022)							(0.028)
time	-0.00001	-0.00002	-0.00003	-0.00002	-0.00001	-0.00001	-0.00001	-0.00001	-0.00002	-0.00003	-0.00002	-0.00001	-0.00001	-0.00001	-0.00003***	-0.00003***	-0.00002**	-0.00001***	-0.00001***	-0.00001***
	(0.00001)	(0.00001)	(0.00002)	(0.00001)	(0.00002)	(0.00001)	(0.00002)	(0.00002)	(0.00002)	(0.00002)	(0.00002)	(0.00002)	(0.00002)	(0.00002)	(0.00001)	(0.00001)	(0.00001)	(0.00000)	(0.00000)	(0.00000)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480
R2	0.002	0.011	0.008	0.006	0.003	0.004	0.008	0.0005	0.013	0.008	0.007	0.004	0.003	0.014	0.074	0.020	0.068	0.715	0.770	0.707
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 1983-1999


mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ cpi + time, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ cpi + time, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ cpi + time, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ cpi + time, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_cpi + time, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ median_cpi + time, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_cpi + time, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ median_cpi + time, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ median_pce + time, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(small ~ median_pce + time, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_pce + time, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ trimmed_cpi + time, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ trimmed_cpi + time, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ trimmed_cpi + time, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(large ~ total_median_shock + time, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(small ~ total_median_shock + time, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ total_median_shock + time, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ total_median_shock + time, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(large ~ pce_def_cycle + time, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_cycle + time, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ pce_def_cycle + time, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(tbills ~ pce_def_cycle + time, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(large ~ pce_def_trend + time, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(small ~ pce_def_trend + time, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_1983_1999_time_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-3.734**							-4.635**							-2.494*			0.093*
	(1.877)							(2.289)							(1.510)			(0.054)
median_cpi		-2.635							-7.854							-0.777			0.491***
		(4.333)							(5.220)							(2.909)			(0.190)
median_pce			-0.181							-7.619
			(11.195)							(12.370)
trimmed_cpi				-2.592							-5.233
				(3.744)							(4.522)
total_median_shock					-4.072**							-4.163					-2.911**			0.019
					(2.071)							(2.538)					(1.267)			(0.047)
pce_def_trend						-0.053**							-0.046
						(0.025)							(0.031)
pce_def_cycle							-0.032							-0.041
							(0.041)							(0.052)
time	-0.00000	0.00000	0.00003	0.00000	0.00003	0.00001	0.00003	-0.00002	-0.0001	-0.0001	-0.00003	0.00002	0.00000	0.00002	-0.0001	-0.00004	-0.00003	-0.00002***	-0.00001***	-0.00002***
	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.0001)	(0.00004)	(0.00004)	(0.00003)	(0.00000)	(0.00000)	(0.00000)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204
R2	0.025	0.003	0.001	0.004	0.024	0.023	0.005	0.025	0.011	0.002	0.007	0.016	0.012	0.005	0.028	0.006	0.031	0.442	0.478	0.432
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 2000 - 2022



mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ cpi + time, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ cpi + time, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ cpi + time, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ cpi + time, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_cpi + time, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ median_cpi + time, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_cpi + time, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(tbills ~ median_cpi + time, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(large ~ median_pce + time, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ median_pce + time, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_pce + time, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(large ~ trimmed_cpi + time, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(small ~ trimmed_cpi + time, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ trimmed_cpi + time, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(large ~ total_median_shock + time, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(small ~ total_median_shock + time, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ total_median_shock + time, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ total_median_shock + time, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_cycle + time, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(small ~ pce_def_cycle + time, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(large ~ pce_def_trend + time, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(small ~ pce_def_trend + time, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_2000_2022_time_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	0.362							1.341							-3.014***			0.058***
	(1.323)							(1.819)							(0.911)			(0.004)
median_cpi		-7.821***							-7.725**							-3.599*			0.580***
		(2.572)							(3.616)							(1.868)			(0.191)
median_pce			-11.523*							-11.212*
			(6.058)							(5.940)
trimmed_cpi				-5.128**							-5.020
				(2.314)							(3.241)
total_median_shock					1.419							2.552					-3.068***			-0.006
					(1.359)							(1.888)					(0.645)			(0.015)
pce_def_trend						-0.005							-0.010
						(0.010)							(0.016)
pce_def_cycle							0.031							0.059*
							(0.024)							(0.032)
time	0.0001	0.0001**	0.0001*	0.0001**	0.0001	0.0001*	0.0001	-0.00001	0.00002	0.00000	0.00002	-0.00000	0.00000	-0.00001	-0.00004	-0.00003	-0.00005**	-0.00001***	-0.00001***	-0.00001***
	(0.00004)	(0.00003)	(0.00004)	(0.00003)	(0.00004)	(0.00003)	(0.00004)	(0.0001)	(0.00005)	(0.00005)	(0.00005)	(0.0001)	(0.0001)	(0.0001)	(0.00002)	(0.00002)	(0.00002)	(0.00000)	(0.00000)	(0.00000)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276
R2	0.011	0.043	0.036	0.028	0.019	0.011	0.026	0.005	0.016	0.013	0.009	0.014	0.002	0.030	0.098	0.026	0.088	0.286	0.435	0.271
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: IV using Ibbotson

# IV regressions


# CPI

mod_large_cpi <- 
  df_common_mo %>% 
ivreg(large ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
ivreg(small ~ cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
ivreg(large ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
ivreg(small ~ median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
ivreg(large ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
ivreg(small ~ median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
ivreg(large ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
ivreg(small ~ trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  ivreg(large ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
  ivreg(small ~ total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
ivreg(large ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
ivreg(small ~ pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
ivreg(large ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
ivreg(small ~ pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-6.239*							-9.206**
	(3.679)							(4.491)
median_cpi		-7.892*							-11.647**
		(4.462)							(5.305)
median_pce			-12.621*							-18.624**
			(7.301)							(8.701)
trimmed_cpi				-6.403*							-9.450**
				(3.630)							(4.312)
total_median_shock					-29.774							-43.938
					(34.556)							(49.035)
pce_def_trend						-0.055*							-0.082**
						(0.033)							(0.040)
pce_def_cycle							0.329							0.485
							(0.345)							(0.472)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	479	479	479	479	479	479	479	479	479	479	479	479	479	479
R2	-0.136	-0.005	-0.027	-0.010	-2.758	-0.052	-1.053	-0.192	0.001	-0.026	-0.006	-3.603	-0.065	-1.259
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

# I have one instrument, so the coefficients are exactly identified (SW p 448) so I can't test for exogeneity statistically. I must assume that the IV is exogenous (that expectations change a year ago (lag(mich_chg, 12) is unrelated to stock returns in the current year, which seems reasonable since stocks are forward looking and their returns today anticipate conditions tomorrow). 
# SEs can be HC only
# To check relevance, examine the F stat of the first stage regression. When there is only one instrument in this regression, the first-stage F-statistic is the square of the t-statistic testing that the coefficient on the instrumental variable. As a rule of thumb, > 10 => instruments not weak so we can rely on the standard methods for statistical inference (SW p 453). In theory, inflation expectations yesterday will influence inflation today under adaptive or hybrid models of inflation expectations formation. This, even if the first stage results are weak, this is a defensible idea. 

# First stage test (use hetero robust errors) for lag(mich_chg, 12)

# IV 1983 to 1999

# CPI

mod_large_cpi <- 
  df_common_mo %>%
  filter(between(year(date), 1983, 1999)) %>% 
ivreg(large ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(small ~ cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(large ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(small ~ median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(large ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(small ~ median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(large ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(small ~ trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
  ivreg(large ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
  ivreg(small ~ total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(large ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(small ~ pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(large ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(small ~ pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_1983_1999_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-3.097							-4.515
	(4.808)							(5.847)
median_cpi		-4.456							-6.497
		(7.090)							(8.843)
median_pce			-4.875							-7.107
			(7.665)							(9.460)
trimmed_cpi				-3.593							-5.238
				(5.703)							(7.132)
total_median_shock					-10.154							-14.803
					(16.223)							(19.275)
pce_def_trend						-0.050							-0.073
						(0.076)							(0.093)
pce_def_cycle							0.586							0.854
							(1.718)							(2.379)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	203	203	203	203	203	203	203	203	203	203	203	203	203	203
R2	0.024	0.003	0.001	0.004	-0.030	0.022	-1.509	0.024	0.011	0.0004	0.007	-0.094	0.007	-2.181
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

# IV 2000-2022


# CPI

mod_large_cpi <- 
  df_common_mo %>%
  filter(between(year(date), 2000, 2022)) %>% 
ivreg(large ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(large ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(large ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>%  
ivreg(large ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
  ivreg(large ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
  ivreg(small ~ total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(large ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(large ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(small ~ pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_2000-2022_Ibbotson.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-8.217							-11.519*
	(5.178)							(6.499)
median_cpi		-10.013*							-14.037**
		(5.704)							(6.858)
median_pce			-19.911*							-27.914*
			(11.865)							(14.704)
trimmed_cpi				-8.102*							-11.358**
				(4.626)							(5.542)
total_median_shock					-45.803							-64.213
					(82.176)							(114.357)
pce_def_trend						-0.059							-0.083*
						(0.037)							(0.044)
pce_def_cycle							0.329							0.461
							(0.366)							(0.483)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	275	275	275	275	275	275	275	275	275	275	275	275	275	275
R2	-0.372	0.011	0.009	-0.004	-9.329	-0.125	-1.425	-0.416	0.004	-0.016	-0.007	-9.669	-0.103	-1.344
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

Regressions: OLS using Russell

# CPI

mod_large_cpi <- 
  df_common_mo %>% 
lm(russell1 ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
lm(russell2 ~  cpi, data = .) 

mod_ltbond_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
lm(tbills ~ cpi, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
lm(russell1 ~ median_cpi, data = .) 

mod_small_median_cpi<- 
 df_common_mo %>% 
lm(russell2 ~  median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
lm(russell1 ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
lm(russell2 ~  median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
lm(russell1 ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
lm(russell2 ~  trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  lm(russell1 ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  lm(russell2 ~  total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
lm(russell1 ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
lm(russell2 ~  pce_def_cycle, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_cycle, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
lm(russell1 ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
lm(russell2 ~  pce_def_trend, data = .) 

mod_ltbond_pce_def_trend <- 
  df_common_mo %>% 
lm(lt_gov ~ pce_def_trend, data = .) 

mod_tbills_pce_def_trend <- 
  df_common_mo %>% 
lm(tbills ~ pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_OLS_russell.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-0.205							-0.406							-2.759***			0.153
	(1.209)							(1.566)							(0.777)			(0.094)
median_cpi		-3.715*							-6.282**							-1.731			1.054*
		(1.952)							(2.477)							(1.356)			(0.549)
median_pce			-2.659							-4.223
			(3.120)							(3.470)
trimmed_cpi				-2.340							-4.455**
				(1.746)							(2.217)
total_median_shock					0.432							0.659					-3.036***			-0.008
					(1.302)							(1.681)					(0.565)			(0.062)
pce_def_trend						-0.010							-0.018
						(0.009)							(0.014)
pce_def_cycle							0.024							0.034
							(0.023)							(0.028)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480	480
R2	0.0002	0.008	0.002	0.004	0.001	0.003	0.006	0.0004	0.013	0.003	0.008	0.001	0.005	0.008	0.057	0.003	0.057	0.030	0.216	0.0001
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 1983-1999


mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell1 ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell2 ~  cpi, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ cpi, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell1 ~ median_cpi, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell2 ~  median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell1 ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(russell2 ~  median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(russell1 ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell2 ~  trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(russell1 ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
  lm(russell2 ~  total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(russell1 ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell2 ~  pce_def_cycle, data = .) 

mod_ltbond_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(lt_gov ~ pce_def_cycle, data = .) 

mod_tbills_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>%  
lm(tbills ~ pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell1 ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 1983, 1999)) %>% 
lm(russell2 ~  pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_1983_1999_OLS_russell.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	-3.996**							-5.060**							-2.035			0.235***
	(1.838)							(2.270)							(1.451)			(0.084)
median_cpi		-3.331							-5.313							0.939			1.036***
		(3.469)							(4.326)							(2.334)			(0.153)
median_pce			-2.366							-1.041
			(4.596)							(5.686)
trimmed_cpi				-3.142							-4.171
				(3.058)							(3.819)
total_median_shock					-4.239**							-5.037*					-2.953**			-0.004
					(2.074)							(2.573)					(1.266)			(0.071)
pce_def_trend						-0.052**							-0.055*
						(0.024)							(0.031)
pce_def_cycle							-0.039							-0.054
							(0.041)							(0.054)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204	204
R2	0.028	0.005	0.001	0.005	0.024	0.022	0.006	0.029	0.007	0.0002	0.006	0.022	0.016	0.007	0.016	0.001	0.026	0.072	0.319	0.00001
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

# OLS 2000 - 2022



mod_large_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(russell1 ~ cpi, data = .) 

mod_small_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(russell2 ~  cpi, data = .) 

mod_ltbond_cpi <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ cpi, data = .) 

mod_tbills_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ cpi, data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(russell1 ~ median_cpi, data = .) 

mod_small_median_cpi<- 
df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(russell2 ~  median_cpi, data = .) 

mod_ltbond_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_cpi, data = .) 

mod_tbills_median_cpi <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(tbills ~ median_cpi, data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(russell1 ~ median_pce, data = .) 

mod_small_median_pce<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(russell2 ~  median_pce, data = .) 

mod_ltbond_median_pce <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ median_pce, data = .) 

# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(russell1 ~ trimmed_cpi, data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(russell2 ~  trimmed_cpi, data = .) 

mod_ltbond_trimmed <- 
  df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(lt_gov ~ trimmed_cpi, data = .) 

# Total median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(russell1 ~ total_median_shock, data = .)

mod_small_median_shock <-  
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
  lm(russell2 ~  total_median_shock, data = .)

mod_ltbond_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>%
lm(lt_gov ~ total_median_shock, data = .) 


mod_tbills_median_shock <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(tbills ~ total_median_shock, data = .) 

# PCE_def_cycle


mod_large_pce_def_cycle <- 
  df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(russell1 ~ pce_def_cycle, data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>% 
lm(russell2 ~  pce_def_cycle, data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
 df_common_mo %>% 
  filter(between(year(date), 2000, 2022)) %>% 
lm(russell1 ~ pce_def_trend, data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
 filter(between(year(date), 2000, 2022)) %>%
lm(russell2 ~  pce_def_trend, data = .) 


rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcov(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcov(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcov(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcov(mod_small_median_cpi))),
              sqrt(diag(vcov(mod_small_median_pce))),
              sqrt(diag(vcov(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_ltbond_cpi ))),
              sqrt(diag(vcov(mod_ltbond_median_cpi))),
              sqrt(diag(vcov(mod_ltbond_median_shock))),
              sqrt(diag(vcovHAC(mod_tbills_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_cpi))),
              sqrt(diag(vcovHAC(mod_tbills_median_shock))))
       

regressions <- 
  list(mod_large_cpi,
  mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock, 
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle,
  mod_ltbond_cpi,
  mod_ltbond_median_cpi,
  mod_ltbond_median_shock,
  mod_tbills_cpi,
  mod_tbills_median_cpi,
  mod_tbills_median_shock
 )

se_list <- 
  list(c("Standard errors", "HC", "IID", "HC", "IID", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "HC", "HC", "HC", "HC", "IID", "IID", "HAC", "HAC", "HAC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_2000_2022_OLS_russell.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_list)

	Dependent variable
	Large cap							Small cap							LT gov. bond			T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)
cpi	0.473							0.592							-3.083***			0.040
	(1.355)							(1.757)							(0.908)			(0.037)
median_cpi		-6.555***							-8.737***							-4.249**			0.351
		(2.541)							(3.289)							(1.805)			(0.287)
median_pce			-11.538*							-12.775**
			(6.181)							(5.598)
trimmed_cpi				-4.077*							-6.112**
				(2.283)							(2.951)
total_median_shock					1.452							1.890					-3.045***			-0.001
					(1.416)							(1.844)					(0.648)			(0.034)
pce_def_trend						-0.003							-0.013
						(0.010)							(0.016)
pce_def_cycle							0.034							0.048
							(0.025)							(0.031)
Standard errors	HC	IID	HC	IID	HC	IID	HC	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HAC	HAC	HAC
Observations	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276	276
R2	0.001	0.024	0.026	0.012	0.009	0.0004	0.019	0.001	0.025	0.019	0.015	0.009	0.003	0.022	0.090	0.020	0.075	0.007	0.064	0.00000
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: IV using Russell

# IV regressions


# CPI

mod_large_cpi <- 
  df_common_mo %>% 
ivreg(russell1 ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
ivreg(russell2 ~  cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
  df_common_mo %>% 
ivreg(russell1 ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
ivreg(russell2 ~  median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
ivreg(russell1 ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
ivreg(russell2 ~  median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
ivreg(russell1 ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
ivreg(russell2 ~  trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
  ivreg(russell1 ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
  ivreg(russell2 ~  total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
ivreg(russell1 ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
ivreg(russell2 ~  pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
ivreg(russell1 ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
ivreg(russell2 ~  pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_russell.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-6.757*							-9.029**
	(3.742)							(4.481)
median_cpi		-8.548*							-11.423**
		(4.507)							(5.344)
median_pce			-13.669*							-18.266**
			(7.396)							(8.790)
trimmed_cpi				-6.936*							-9.268**
				(3.666)							(4.346)
total_median_shock					-32.248							-43.093
					(36.883)							(48.038)
pce_def_trend						-0.060*							-0.080**
						(0.034)							(0.040)
pce_def_cycle							0.356							0.476
							(0.365)							(0.469)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	479	479	479	479	479	479	479	479	479	479	479	479	479	479
R2	-0.155	-0.004	-0.030	-0.010	-3.171	-0.060	-1.219	-0.166	0.006	-0.027	-0.0002	-3.521	-0.055	-1.336
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

# I have one instrument, so the coefficients are exactly identified (SW p 448) so I can't test for exogeneity statistically. I must assume that the IV is exogenous (that expectations change a year ago (lag(mich_chg, 12) is unrelated to stock returns in the current year, which seems reasonable since stocks are forward looking and their returns today anticipate conditions tomorrow). 
# SEs can be HC only
# To check relevance, examine the F stat of the first stage regression. When there is only one instrument in this regression, the first-stage F-statistic is the square of the t-statistic testing that the coefficient on the instrumental variable. As a rule of thumb, > 10 => instruments not weak so we can rely on the standard methods for statistical inference (SW p 453). In theory, inflation expectations yesterday will influence inflation today under adaptive or hybrid models of inflation expectations formation. This, even if the first stage results are weak, this is a defensible idea. 

# First stage test (use hetero robust errors) for lag(mich_chg, 12)

# IV 1983 to 1999

# CPI

mod_large_cpi <- 
  df_common_mo %>%
  filter(between(year(date), 1983, 1999)) %>% 
ivreg(russell1 ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(russell2 ~  cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell1 ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell2 ~  median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(russell1 ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell2 ~  median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
ivreg(russell1 ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell2 ~  trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
  ivreg(russell1 ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>% 
  ivreg(russell2 ~  total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell1 ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell2 ~  pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell1 ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
filter(between(year(date), 1983, 1999)) %>%
ivreg(russell2 ~  pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_1983_1999.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV_russell", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-2.847							-2.344
	(4.808)							(6.333)
median_cpi		-4.097							-3.373
		(7.106)							(9.403)
median_pce			-4.482							-3.690
			(7.684)							(10.188)
trimmed_cpi				-3.303							-2.720
				(5.716)							(7.576)
total_median_shock					-9.335							-7.686
					(15.929)							(19.807)
pce_def_trend						-0.046							-0.038
						(0.076)							(0.102)
pce_def_cycle							0.539							0.443
							(1.649)							(1.763)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	203	203	203	203	203	203	203	203	203	203	203	203	203	203
R2	0.026	0.005	0.001	0.006	-0.012	0.021	-1.314	0.021	0.007	0.0003	0.006	0.013	0.013	-0.631
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

# IV 2000-2022


# CPI

mod_large_cpi <- 
  df_common_mo %>%
  filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell1 ~ cpi | lag(mich), data = .) 


mod_small_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  cpi | lag(mich), data = .) 


# Median CPI

mod_large_median_cpi <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell1 ~ median_cpi | lag(mich), data = .) 


mod_small_median_cpi<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  median_cpi | lag(mich), data = .) 

# Median PCE

mod_large_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell1 ~ median_pce | lag(mich), data = .) 

mod_small_median_pce <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  median_pce | lag(mich), data = .) 


# Trimmed mean CPI

mod_large_trimmed <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>%  
ivreg(russell1 ~ trimmed_cpi | lag(mich), data = .) 

mod_small_trimmed<- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  trimmed_cpi | lag(mich), data = .) 

# Median shock

mod_large_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
  ivreg(russell1 ~ total_median_shock | lag(mich), data = .)

mod_small_median_shock <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
  ivreg(russell2 ~  total_median_shock | lag(mich), data = .)


# PCE_def_cycle

mod_large_pce_def_cycle <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell1 ~ pce_def_cycle | lag(mich), data = .) 

mod_small_pce_def_cycle <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  pce_def_cycle | lag(mich), data = .) 

# PCE_def_trend


mod_large_pce_def_trend <- 
  df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell1 ~ pce_def_trend | lag(mich), data = .) 

mod_small_pce_def_trend <- 
 df_common_mo %>% 
filter(between(year(date), 2000, 2022)) %>% 
ivreg(russell2 ~  pce_def_trend | lag(mich), data = .) 


rob_se <- 
  list(       sqrt(diag(vcovHC(mod_large_cpi))),
              sqrt(diag(vcovHC(mod_large_median_cpi))),
              sqrt(diag(vcovHC(mod_large_median_pce))),
              sqrt(diag(vcovHC(mod_large_trimmed))),
              sqrt(diag(vcovHC(mod_large_median_shock))),
              sqrt(diag(vcovHC(mod_large_pce_def_trend))),
              sqrt(diag(vcovHC(mod_large_pce_def_cycle))),
              sqrt(diag(vcovHC(mod_small_cpi))),
              sqrt(diag(vcovHC(mod_small_median_cpi))),
              sqrt(diag(vcovHC(mod_small_median_pce))),
              sqrt(diag(vcovHC(mod_small_trimmed))),
              sqrt(diag(vcovHC(mod_small_median_shock))),
              sqrt(diag(vcovHC(mod_small_pce_def_trend))),
              sqrt(diag(vcovHC(mod_small_pce_def_cycle))))
       

regressions <- 
  list(mod_large_cpi,
       mod_large_median_cpi,
  mod_large_median_pce,
  mod_large_trimmed,
  mod_large_median_shock,
  mod_large_pce_def_trend,
  mod_large_pce_def_cycle,
  mod_small_cpi,
  mod_small_median_cpi,
  mod_small_median_pce,
  mod_small_trimmed,
  mod_small_median_shock,
  mod_small_pce_def_trend,
  mod_small_pce_def_cycle)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "main_regressions_IV_2000-2022_russell.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted", "Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations"), 
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap IV", "Small cap IV"),
  model.numbers = T,
  add.lines = list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC")))

	Dependent variable
	Large cap IV							Small cap IV
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
cpi	-9.036*							-12.246*
	(5.342)							(6.558)
median_cpi		-11.011*							-14.923**
		(5.777)							(6.806)
median_pce			-21.896*							-29.676**
			(12.112)							(14.679)
trimmed_cpi				-8.909*							-12.075**
				(4.681)							(5.504)
total_median_shock					-50.370							-68.267
					(89.803)							(120.901)
pce_def_trend						-0.065*							-0.088**
						(0.038)							(0.044)
pce_def_cycle							0.361							0.490
							(0.392)							(0.508)
Standard errors	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC	HC
Observations	275	275	275	275	275	275	275	275	275	275	275	275	275	275
R2	-0.432	0.011	0.005	-0.006	-10.885	-0.141	-1.679	-0.466	0.012	-0.014	0.0004	-11.851	-0.121	-1.817
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted
	Instrument is lag of University of Michigan Year-Ahead Consumer Survey Inflation Expectations

Figure 4: rolling inflation betas

roll_df <- 
  df_common_mo %>% 
  mutate(lag_median_cpi = median_cpi - lag(median_cpi),
         lag_total_shock = total_median_shock - lag(total_median_shock))


plot1 <- 
  roll_lm(x = roll_df$median_cpi, y = roll_df$lag_median_cpi, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1984-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  annotate("text", x = as.Date("2008-01-01"), y = 1.3, label = "Coefficient =\n 0.73") +
  theme +
  labs(x = NULL, y = "Median CPI autoregression coefficient") +
  scale_x_date(expand = c(0, 0))


plot2 <- 
  roll_lm(x = roll_df$total_median_shock, y = roll_df$lag_total_shock, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1984-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  annotate("text", x = as.Date("2008-01-01"), y = 1.1, label = "Coefficient =\n 0.42") +
  theme +
  labs(x = NULL, y = "Total median CPI shock autoregression coefficient") +
  scale_x_date(expand = c(0, 0))

plot1 + plot2 + plot_annotation(caption = "Source: FRED")

# for FO AC

df_common_mo %>% 
  lm(total_median_shock ~ lag(total_median_shock, 1), data = .) %>% summary

## 
## Call:
## lm(formula = total_median_shock ~ lag(total_median_shock, 1), 
##     data = .)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0148033 -0.0009606  0.0000223  0.0010419  0.0087211 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                -8.132e-05  1.001e-04  -0.812    0.417    
## lag(total_median_shock, 1)  4.244e-01  4.159e-02  10.205   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002187 on 477 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.1792, Adjusted R-squared:  0.1775 
## F-statistic: 104.1 on 1 and 477 DF,  p-value: < 2.2e-16

# Rolling regressions


plot1 <- 
  roll_lm(x = df_common_mo$median_cpi, y = df_common_mo$large, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1988-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  theme +
  labs(x = NULL, y = "Large cap median") +
  scale_x_date(expand = c(0, 0))

plot2 <- 
  roll_lm(x = df_common_mo$total_median_shock, y = df_common_mo$large, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1988-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  theme +
  labs(x = NULL, y = "Large cap headline") +
  scale_x_date(expand = c(0, 0))
  
  # <- median CPI

plot3 <- 
  roll_lm(x = df_common_mo$median_cpi, y = df_common_mo$small, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1988-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  theme +
  labs(x = NULL, y = "Small cap median") +
  scale_x_date(expand = c(0, 0))


plot4 <- 
  roll_lm(x = df_common_mo$total_median_shock, y = df_common_mo$small, width = 60, min_obs = 60) %>% 
  as_tibble %>% 
  select(1,3) %>% 
  set_colnames(c("one", "three")) %>% 
  mutate(beta = one[,2], se = three[,2]) %>% 
  drop_na %>% 
  mutate(date = seq(as.Date("1988-01-01"),
                    by = "month",
                    length = length(beta))) %>% 
  select(date, beta, se) %>% 
  ggplot(aes(x = date, y = beta)) +
  geom_ribbon(aes(ymin = beta - 2*se, ymax = beta + 2*se), fill = "lightgray") +
  geom_line() +
  theme +
  labs(x = NULL, y = "Small cap headline") +
  scale_x_date(expand = c(0, 0))

(plot1 + plot2) / (plot3 + plot4) + plot_annotation(caption = "Source: CFA Institute, FRED") 

Regressions: Phillips curve (economic activity and inflation expectations)

# payrolls

payrolls_df <- 
sbbi_mo %>% 
  left_join(payrolls, by = "date") %>% 
  left_join(tightening, by = "date") 


mod_large_payrolls <- 
  payrolls_df %>% 
  lm(large ~ payrolls_cycle + tightening, data = .)

mod_small_payrolls <- 
  payrolls_df %>% 
  lm(small ~ payrolls_cycle + tightening, data = .)

mod_ltbond_payrolls <- 
  payrolls_df %>% 
  lm(lt_gov ~ payrolls_cycle + tightening, data = .)

mod_tbills_payrolls <- 
  payrolls_df %>% 
  lm(tbills ~ payrolls_cycle + tightening, data = .)


# rpce

rpce_df <- 
sbbi_mo %>% 
  right_join(rpce, by = "date") %>% 
  left_join(tightening, by = "date")

mod_large_rpce <- 
  rpce_df %>% 
  lm(large ~ rpce_cycle + tightening, data = .)

mod_small_rpce <- 
  rpce_df %>% 
  lm(small ~ rpce_cycle + tightening, data = .)

mod_ltbond_rpce <- 
  rpce_df %>% 
  lm(lt_gov ~ rpce_cycle + tightening, data = .)

mod_tbills_rpce <- 
  rpce_df %>% 
  lm(tbills ~ rpce_cycle + tightening, data = .)

# hot tcu


tcu_df <- 
sbbi_mo %>% 
  right_join(tcu, by = "date") %>% 
  right_join(tightening, by = "date")

mod_large_tcu <- 
  tcu_df %>% 
  lm(large ~ hot_tcu + tightening, data = .)

mod_small_tcu <- 
  tcu_df %>% 
  lm(small~ hot_tcu + tightening, data = .)

mod_ltbond_tcu <- 
  tcu_df %>% 
  lm(lt_gov ~ hot_tcu + tightening, data = .)

mod_tbills_tcu <- 
  tcu_df %>% 
  lm(tbills ~ hot_tcu + tightening, data = .)

# expectations df

exp_infl <- 
c("MICH", "EXPINF1YR","CPIAUCSL", "PCEPI") %>% 
  tq_get(get = "economic.data", from = "1947-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  drop_na %>% 
  select(date, mich = MICH, cleve = EXPINF1YR, cpi = CPIAUCSL, pce = PCEPI) %>% 
  mutate(across(2:3, ~ ((1 + (./100))^(1/12))) - 1, #convert mich and cleve to monthly 
         across(4:5, ~ (./lag(., 1)) - 1)) 

exp_infl_df <- 
  exp_infl %>% 
  mutate(mich_cpi_surprise = mich - cpi, 
         cleve_cpi_surprise = cleve - cpi,
         mich_pce_surprise = mich - pce,
         cleve_pce_surprise = cleve - pce) %>% 
  left_join(sbbi_mo, by = "date") %>% 
  left_join(be_infl, by = "date") %>%   
  mutate(be_infl = ((1 + (be_infl/100))^(1/12)) - 1) 

# mich

mod_large_mich <- 
  exp_infl_df %>% 
  lm(large ~ mich, data = .)

mod_small_mich <- 
  exp_infl_df %>% 
  lm(small ~ mich, data = .)

mod_ltbond_mich <- 
  exp_infl_df %>% 
  lm(lt_gov ~ mich, data = .)

mod_tbills_mich <- 
  exp_infl_df %>% 
  lm(tbills ~ mich, data = .)

#cleve
  
mod_large_cleve <- 
  exp_infl_df %>% 
  lm(large ~ cleve, data = .)

mod_small_cleve <- 
  exp_infl_df %>% 
  lm(small ~ cleve, data = .)

mod_ltbond_cleve <- 
  exp_infl_df %>% 
  lm(lt_gov ~ cleve, data = .)

mod_tbills_cleve <- 
  exp_infl_df %>% 
  lm(tbills ~ cleve, data = .)

#be infl
  
mod_large_be <- 
  exp_infl_df %>% 
  lm(large ~ be_infl, data = .)

mod_small_be <- 
  exp_infl_df %>% 
  lm(small ~ be_infl, data = .)

mod_ltbond_be <- 
  exp_infl_df %>% 
  lm(lt_gov ~ be_infl, data = .)

mod_tbills_be <- 
  exp_infl_df %>% 
  lm(tbills ~ be_infl, data = .)


# mich sur

mod_large_mich_cpi_sur <- 
  exp_infl_df %>% 
  lm(large ~ mich_cpi_surprise, data = .)

mod_small_mich_cpi_sur <- 
  exp_infl_df %>% 
  lm(small ~ mich_cpi_surprise, data = .)

mod_ltbond_mich_cpi_sur <- 
  exp_infl_df %>% 
  lm(lt_gov ~ mich_cpi_surprise, data = .)

mod_tbills_mich_cpi_sur <- 
  exp_infl_df %>% 
  lm(tbills ~ mich_cpi_surprise, data = .)
    

    # regressions



rob_se <- 
  list(sqrt(diag(vcov(mod_large_payrolls))),
       sqrt(diag(vcovHC(mod_large_rpce))),
       sqrt(diag(vcov(mod_large_tcu))),
       sqrt(diag(vcov(mod_large_mich))),
       sqrt(diag(vcov(mod_large_cleve))),
       sqrt(diag(vcovHC(mod_large_mich_cpi_sur))),
       sqrt(diag(vcovHC(mod_small_payrolls))),
       sqrt(diag(vcovHC(mod_small_rpce))),
       sqrt(diag(vcovHC(mod_small_tcu))),
       sqrt(diag(vcov(mod_small_mich))),
       sqrt(diag(vcov(mod_small_cleve))),
       sqrt(diag(vcovHC(mod_small_mich_cpi_sur))),
       sqrt(diag(vcovHC(mod_ltbond_payrolls))),
       sqrt(diag(vcovHC(mod_ltbond_rpce))),
       sqrt(diag(vcovHC(mod_ltbond_tcu))),
       sqrt(diag(vcov(mod_ltbond_mich))),
       sqrt(diag(vcovHC(mod_ltbond_cleve))),
       sqrt(diag(vcovHC(mod_ltbond_mich_cpi_sur))),
       sqrt(diag(vcovHAC(mod_tbills_payrolls))),
       sqrt(diag(vcovHAC(mod_tbills_rpce))),
       sqrt(diag(vcovHAC(mod_tbills_tcu))),
       sqrt(diag(vcovHAC(mod_tbills_mich))),
       sqrt(diag(vcovHAC(mod_tbills_cleve))),
       sqrt(diag(vcovHAC(mod_tbills_mich_cpi_sur))))
       
    se_notes <- 
  list(c("Standard errors", "IID", "HC", "IID", "IID", "IID", "HC","HC", "HC", "HC", "IID", "IID",  "HC", "HC", "HC", "HC", "IID", "HC", "HC", "IID", "IID", "IID", "IID", "IID", "IID"))
       

regressions <- 
  list(mod_large_payrolls,
       mod_large_rpce,
       mod_large_tcu,
       mod_large_mich,
       mod_large_cleve,
       mod_large_mich_cpi_sur,
       mod_small_payrolls,
       mod_small_rpce,
       mod_small_tcu,
       mod_small_mich,
       mod_small_cleve,
       mod_small_mich_cpi_sur,
       mod_ltbond_payrolls,
       mod_ltbond_rpce,
       mod_ltbond_tcu,
       mod_ltbond_mich,
       mod_ltbond_cleve,
       mod_ltbond_mich_cpi_sur,
       mod_tbills_payrolls,
       mod_tbills_rpce,
       mod_tbills_tcu,
       mod_tbills_mich,
       mod_tbills_cleve,
       mod_tbills_mich_cpi_sur)

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq" ),
  omit = c("Constant"),
  type = "html",
  out = "Phillips_regressions.html",
  notes = c("Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_notes)

	Dependent variable
	Large cap						Small cap						LT gov. bond						T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)	(18)	(19)	(20)	(21)	(22)	(23)	(24)
payrolls_cycle	-0.00000***						-0.00001***						0.00000						0.00000*
	(0.00000)						(0.00000)						(0.00000)						(0.00000)
rpce_cycle		-0.0001						-0.0001						-0.00003						-0.00000
		(0.00005)						(0.0001)						(0.00002)						(0.00001)
hot_tcu1			-0.002						-0.004						-0.004						0.001**
			(0.004)						(0.005)						(0.003)						(0.001)
tightening1	-0.004	-0.007*	-0.008**				-0.004	-0.008*	-0.008				-0.003	-0.002	-0.002				0.001**	0.002***	0.002***
	(0.003)	(0.003)	(0.004)				(0.004)	(0.005)	(0.005)				(0.002)	(0.002)	(0.003)				(0.001)	(0.001)	(0.001)
mich				-5.674						-10.346**						-1.574						1.116
				(3.509)						(4.555)						(2.488)						(0.756)
cleve					-0.566						-2.373						3.278*						2.265***
					(2.039)						(2.652)						(1.712)						(0.096)
mich_cpi_surprise						0.157						-0.575						3.056***						-0.143**
						(1.229)						(1.631)						(0.798)						(0.061)
Standard errors	IID	HC	IID	IID	IID	HC	HC	HC	HC	IID	IID	HC	HC	HC	HC	IID	HC	HC	IID	IID	IID	IID	IID	IID
Observations	828	779	684	503	503	502	828	779	684	503	503	502	828	779	684	503	503	502	828	779	684	503	503	502
Adjusted R2	0.027	0.011	0.005	0.003	-0.002	-0.002	0.031	0.006	0.003	0.008	-0.0004	-0.001	0.0001	0.003	0.003	-0.001	0.008	0.059	0.129	0.105	0.185	0.062	0.784	0.019
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Figure 5: food and energy shocks

# Figure __. Line chart of food and energy shocks

cpi_short <- 
  c("CPIAUCSL", "CPILFESL") %>% 
  tq_get(get = "economic.data", from = "1957-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  mutate(across(2:3, ~ 100*((./lag(., 12)) - 1))) %>% 
  drop_na

cpi_short %>%
  filter(year(date) < 2023) %>% 
  rename(cpi = cpiaucsl, cpi_core = cpilfesl) %>% 
  mutate(fe_shock = cpi - cpi_core) %>% 
  ggplot(aes(x = date, y = fe_shock)) +
  geom_rect(aes(xmin = as.Date("1970-01-01"), xmax = as.Date("1980-12-01"),
                ymin = -Inf, ymax = +Inf), fill = "lightgray") +
  geom_rect(aes(xmin = as.Date("2005-01-01"), xmax = as.Date("2008-12-01"),
                ymin = -Inf, ymax = +Inf), fill = "lightgray") +
  geom_rect(aes(xmin = as.Date("2021-01-01"), xmax = as.Date("2022-12-01"),
                ymin = -Inf, ymax = +Inf), fill = "lightgray") +
  geom_line() +
  theme +
  labs(x = NULL, y = "12-month headline CPI less 12-month core CPI") +
  scale_x_date(expand = c(0, 0))

Regressions: supply shocks using Ibbotson

# FE shock

mod_large_fe_shock <- 
  df_common_mo %>% 
  lm(large ~ fe_shock + tightening + mich, data = .)

mod_small_fe_shock <- 
  df_common_mo %>% 
  lm(small ~ fe_shock + tightening + mich, data = .)

mod_ltbond_fe_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ fe_shock + tightening + mich, data = .) 

mod_tbills_fe_shock <- 
  df_common_mo %>% 
lm(tbills ~ fe_shock + tightening + mich, data = .) 

# Other shock

mod_large_other_median_shock <- 
  df_common_mo %>% 
  lm(large ~ other_median_shock + tightening + mich, data = .)

mod_small_other_median_shock <- 
  df_common_mo %>% 
  lm(small ~ other_median_shock + tightening + mich, data = .)

mod_ltbond_other_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ other_median_shock + tightening + mich, data = .) 

mod_tbills_other_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ other_median_shock + tightening + mich, data = .) 

# Orders shock

mod_large_ismmsdi <- 
  df_common_mo %>% 
lm(large ~ ismmsdi_diff + tightening + mich, data = .) 

mod_small_ismmsdi <- 
  df_common_mo %>% 
lm(small ~ ismmsdi_diff + tightening + mich, data = .) 

mod_ltbond_ismmsdi <- 
  df_common_mo %>% 
lm(lt_gov ~ ismmsdi_diff + tightening + mich, data = .) 

mod_tbills_ismmsdi <- 
  df_common_mo %>% 
lm(tbills ~ ismmsdi_diff + tightening + mich, data = .) 


# Supply shock

mod_large_supply <- 
  df_common_mo %>% 
lm(large ~ gscpi_diff + tightening + mich, data = .) 

mod_small_supply <- 
  df_common_mo %>% 
lm(small ~ gscpi_diff + tightening + mich, data = .) 

mod_ltbond_supply <- 
  df_common_mo %>% 
lm(lt_gov ~ gscpi_diff + tightening + mich, data = .) 

mod_tbills_supply <- 
  df_common_mo %>% 
lm(tbills ~ gscpi_diff + tightening + mich, data = .) 




rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_fe_shock))),
       sqrt(diag(vcovHC(mod_large_other_median_shock))),
       sqrt(diag(vcov(mod_large_supply))),
       sqrt(diag(vcov(mod_large_ismmsdi))),
       sqrt(diag(vcovHC(mod_small_fe_shock))),
       sqrt(diag(vcovHC(mod_small_other_median_shock))),
       sqrt(diag(vcov(mod_small_ismmsdi))),
       sqrt(diag(vcov(mod_small_supply))),
       sqrt(diag(vcovHC(mod_ltbond_fe_shock))),
       sqrt(diag(vcov(mod_ltbond_other_median_shock))),
       sqrt(diag(vcov(mod_ltbond_ismmsdi))),
       sqrt(diag(vcov(mod_ltbond_supply))),
       sqrt(diag(vcovHAC(mod_tbills_fe_shock))),
       sqrt(diag(vcovHAC(mod_tbills_other_median_shock))),
       sqrt(diag(vcovHAC(mod_tbills_ismmsdi))),
       sqrt(diag(vcovHAC(mod_tbills_supply))))
       
    se_notes <- 
  list(c("Standard errors", "HC", "HC", "IID", "IID", "HC", "HC", "IID", "IID",  "HC", "IID", "IID", "IID",  "HAC", "HAC", "HAC", "HAC"))   
       

regressions <- 
  list(mod_large_fe_shock,
       mod_large_other_median_shock,
       mod_large_ismmsdi,
       mod_large_supply,
       mod_small_fe_shock,
       mod_small_other_median_shock,
       mod_small_ismmsdi,
       mod_small_supply,
       mod_ltbond_fe_shock,
       mod_ltbond_other_median_shock,
       mod_ltbond_ismmsdi,
       mod_ltbond_supply,
       mod_tbills_fe_shock,
       mod_tbills_other_median_shock,
       mod_tbills_ismmsdi,
       mod_tbills_supply)

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq" ),
  omit = c("Constant"),
  type = "html",
  out = "supply_shock_regressions_Ibbotson.html",
  notes = c("Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_notes)

	Dependent variable
	Large cap				Small cap				LT gov. bond				T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)
fe_shock	0.748				2.456				-3.353***				-0.103
	(1.215)				(1.603)				(0.936)				(0.086)
other_median_shock		1.209				-0.429				1.076				0.345
		(2.961)				(3.580)				(1.424)				(0.223)
ismmsdi_diff			-0.001				0.0001				-0.001				-0.00000
							(0.001)				(0.001)				(0.00002)
gscpi_diff				0.014				0.020**				-0.001				-0.00000
								(0.010)				(0.005)				(0.0001)
tightening1	-0.001	-0.001	-0.001	0.001	0.001	0.001	0.001	0.007	0.001	0.001	0.001	0.002	0.001**	0.002**	0.002***	0.002***
	(0.005)	(0.005)	(0.007)	(0.005)	(0.007)	(0.007)	(0.006)	(0.009)	(0.003)	(0.003)	(0.003)	(0.005)	(0.001)	(0.001)	(0.001)	(0.0004)
mich	-0.007*	-0.007*	-0.007*	-0.006**	-0.013***	-0.011**	-0.011***	-0.008*	-0.001	-0.004*	-0.004*	-0.005**	0.0004	0.0003	0.0003	-0.0003
	(0.004)	(0.004)	(0.004)	(0.003)	(0.005)	(0.004)	(0.004)	(0.005)	(0.003)	(0.002)	(0.002)	(0.003)	(0.001)	(0.001)	(0.001)	(0.0003)
Standard errors	HC	HC	IID	IID	HC	HC	IID	IID	HC	IID	IID	IID	HAC	HAC	HAC	HAC
Observations	480	480	480	299	480	480	480	299	480	480	480	299	480	480	480	299
Adjusted R2	0.005	0.004	0.005	0.015	0.018	0.008	0.008	0.016	0.061	0.001	0.004	0.004	0.097	0.108	0.087	0.164
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: supply shocks using Russell

# FE shock

mod_large_fe_shock <- 
  df_common_mo %>% 
  lm(russell1 ~ fe_shock + tightening + mich, data = .)

mod_small_fe_shock <- 
  df_common_mo %>% 
  lm(russell2 ~ fe_shock + tightening + mich, data = .)

mod_ltbond_fe_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ fe_shock + tightening + mich, data = .) 

mod_tbills_fe_shock <- 
  df_common_mo %>% 
lm(tbills ~ fe_shock + tightening + mich, data = .) 

# Other shock

mod_large_other_median_shock <- 
  df_common_mo %>% 
  lm(russell1 ~other_median_shock + tightening + mich, data = .)

mod_small_other_median_shock <- 
  df_common_mo %>% 
  lm(russell2 ~other_median_shock + tightening + mich, data = .)

mod_ltbond_other_median_shock <- 
  df_common_mo %>% 
lm(lt_gov ~ other_median_shock + tightening + mich, data = .) 

mod_tbills_other_median_shock <- 
  df_common_mo %>% 
lm(tbills ~ other_median_shock + tightening + mich, data = .) 

# Orders shock

mod_large_ismmsdi <- 
  df_common_mo %>% 
lm(russell1 ~ismmsdi_diff + tightening + mich, data = .) 

mod_small_ismmsdi <- 
  df_common_mo %>% 
lm(russell2 ~ismmsdi_diff + tightening + mich, data = .) 

mod_ltbond_ismmsdi <- 
  df_common_mo %>% 
lm(lt_gov ~ ismmsdi_diff + tightening + mich, data = .) 

mod_tbills_ismmsdi <- 
  df_common_mo %>% 
lm(tbills ~ ismmsdi_diff + tightening + mich, data = .) 


# Supply shock

mod_large_supply <- 
  df_common_mo %>% 
lm(russell1 ~gscpi_diff + tightening + mich, data = .) 

mod_small_supply <- 
  df_common_mo %>% 
lm(russell2 ~gscpi_diff + tightening + mich, data = .) 

mod_ltbond_supply <- 
  df_common_mo %>% 
lm(lt_gov ~ gscpi_diff + tightening + mich, data = .) 

mod_tbills_supply <- 
  df_common_mo %>% 
lm(tbills ~ gscpi_diff + tightening + mich, data = .) 




rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_fe_shock))),
       sqrt(diag(vcovHC(mod_large_other_median_shock))),
       sqrt(diag(vcov(mod_large_supply))),
       sqrt(diag(vcov(mod_large_ismmsdi))),
       sqrt(diag(vcovHC(mod_small_fe_shock))),
       sqrt(diag(vcovHC(mod_small_other_median_shock))),
       sqrt(diag(vcov(mod_small_ismmsdi))),
       sqrt(diag(vcov(mod_small_supply))),
       sqrt(diag(vcovHC(mod_ltbond_fe_shock))),
       sqrt(diag(vcov(mod_ltbond_other_median_shock))),
       sqrt(diag(vcov(mod_ltbond_ismmsdi))),
       sqrt(diag(vcov(mod_ltbond_supply))),
       sqrt(diag(vcovHAC(mod_tbills_fe_shock))),
       sqrt(diag(vcovHAC(mod_tbills_other_median_shock))),
       sqrt(diag(vcovHAC(mod_tbills_ismmsdi))),
       sqrt(diag(vcovHAC(mod_tbills_supply))))
       
    se_notes <- 
  list(c("Standard errors", "HC", "HC", "IID", "IID", "HC", "HC", "IID", "IID",  "HC", "IID", "IID", "IID",  "HAC", "HAC", "HAC", "HAC"))   
       

regressions <- 
  list(mod_large_fe_shock,
       mod_large_other_median_shock,
       mod_large_ismmsdi,
       mod_large_supply,
       mod_small_fe_shock,
       mod_small_other_median_shock,
       mod_small_ismmsdi,
       mod_small_supply,
       mod_ltbond_fe_shock,
       mod_ltbond_other_median_shock,
       mod_ltbond_ismmsdi,
       mod_ltbond_supply,
       mod_tbills_fe_shock,
       mod_tbills_other_median_shock,
       mod_tbills_ismmsdi,
       mod_tbills_supply)

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "rsq" ),
  omit = c("Constant"),
  type = "html",
  out = "supply_shock_regressions_Russell.html",
  notes = c("Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T,
  add.lines = se_notes)

	Dependent variable
	Large cap				Small cap				LT gov. bond				T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)
fe_shock	0.826				1.674				-3.353***				-0.103
	(1.238)				(1.581)				(0.936)				(0.086)
other_median_shock		1.117				-0.574				1.076				0.345
		(3.031)				(3.562)				(1.424)				(0.223)
ismmsdi_diff			-0.001				-0.00004				-0.001				-0.00000
							(0.001)				(0.001)				(0.00002)
gscpi_diff				0.015				0.020**				-0.001				-0.00000
								(0.010)				(0.005)				(0.0001)
tightening1	-0.0002	-0.0003	-0.001	0.002	-0.0005	-0.001	-0.001	0.006	0.001	0.001	0.001	0.002	0.001**	0.002**	0.002***	0.002***
	(0.005)	(0.005)	(0.007)	(0.005)	(0.006)	(0.006)	(0.006)	(0.009)	(0.003)	(0.003)	(0.003)	(0.005)	(0.001)	(0.001)	(0.001)	(0.0004)
mich	-0.008*	-0.007*	-0.007**	-0.007**	-0.012**	-0.011**	-0.011***	-0.010**	-0.001	-0.004*	-0.004*	-0.005**	0.0004	0.0003	0.0003	-0.0003
	(0.004)	(0.004)	(0.004)	(0.003)	(0.005)	(0.004)	(0.004)	(0.005)	(0.003)	(0.002)	(0.002)	(0.003)	(0.001)	(0.001)	(0.001)	(0.0003)
Standard errors	HC	HC	IID	IID	HC	HC	IID	IID	HC	IID	IID	IID	HAC	HAC	HAC	HAC
Observations	480	480	480	299	480	480	480	299	480	480	480	299	480	480	480	299
Adjusted R2	0.007	0.005	0.006	0.018	0.014	0.009	0.009	0.021	0.061	0.001	0.004	0.004	0.097	0.108	0.087	0.164
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted

Regressions: money supply using Ibbotson

money_cpi <- 
c("M2SL", "CPIAUCSL") %>% 
  tq_get(get = "economic.data", from = "1947-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  filter(year(date) > 1982)

money_cpi %>% 
  mutate(across(2:3, ~ 100*(((./lag(., 12)) - 1)))) %>% 
  #cor.test(~ M2SL + CPIAUCSL, data = .)
  ggplot(aes(x = M2SL, y = CPIAUCSL)) +
  geom_point() +
  theme +
  labs(x = "M2 percentage change", y = "CPI price level change", caption = "Source: FRED") +
  # annotate("text", x = 60, y = 90, label = "Corr. = 0.75", family = "serif", 
  #          size = 8) + ylim(0, 125)
  annotate("text", x = 19, y = 10, label = "Corr. = 0.06", family = "serif",
           size = 8) + ylim(0, 17)

returns <- 
read_csv("sbbi_monthly.csv", locale = locale(encoding = "latin1")) %>% 
    select(c(1, 2, 5, 7, 15)) %>% 
    set_colnames(c("date", "large", "small", "lt_gov", "tbills")) %>% 
  mutate(date = seq(as.Date("1926-01-01"),
                    by = "month",
                    length = length(large))) %>% 
  mutate(across(2:5, ~ cumprod(1 + .))) 


# large

mod_large_short <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 12)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(large ~ m2sl, data = .)


mod_large_long <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 120)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(large ~ m2sl, data = .)

# small

mod_small_short <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 12)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(small ~m2sl, data = .)


mod_small_long <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 120)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(small ~ m2sl, data = .)

# bonds

mod_bonds_short <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 12)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(lt_gov ~ m2sl, data = .)


mod_bonds_long <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 120)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(lt_gov ~ m2sl, data = .)

# tbills

mod_tbills_short <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 12)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(tbills ~ m2sl, data = .)

mod_tbills_long <- 
  returns %>% 
  right_join(money_cpi, by = "date") %>% 
  drop_na %>% 
  mutate(across(2:7, ~ (. / lag(., 120)) - 1)) %>% 
  drop_na %>% 
  clean_names %>% 
  lm(tbills ~ m2sl, data = .)


  rob_se <- 
  list(sqrt(diag(vcovHAC(mod_large_short))),
              sqrt(diag(vcovHAC(mod_large_long))),
              sqrt(diag(vcovHAC(mod_small_short))),
              sqrt(diag(vcovHAC(mod_small_long))),
              sqrt(diag(vcovHAC(mod_bonds_short))),
              sqrt(diag(vcovHAC(mod_bonds_long))),
              sqrt(diag(vcovHAC(mod_tbills_short))),
              sqrt(diag(vcovHAC(mod_tbills_long))))
       

regressions <- 
  list(mod_large_short,
       mod_large_long,
       mod_small_short,
       mod_small_long,
       mod_bonds_short,
       mod_bonds_long,
       mod_tbills_short,
       mod_tbills_long)


stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "money_supply_regressions.html",
  notes = c("Standard errors are heteroskedasticity and autocorrelation adjusted",
            "Source: author's regressions"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap", "LT gov. bond", "T-bills"),
  model.numbers = T)

	Dependent variable
	Large cap		Small cap		LT gov. bond		T-bills
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
m2sl	0.200	-2.950	-0.215	-1.569*	0.470	-1.911***	-0.084*	-1.009***
	(0.582)	(2.301)	(0.792)	(0.811)	(0.333)	(0.222)	(0.045)	(0.076)
Observations	474	366	474	366	474	366	474	366
R2	0.002	0.218	0.002	0.142	0.024	0.483	0.014	0.602
Note:	p<0.1; p<0.05; p<0.01
	Standard errors are heteroskedasticity and autocorrelation adjusted
	Source: author’s regressions

Regressions: inflation vol

# Dataframes from monthly data

cpi_vol_mo <- 
cpi_short_mo %>% 
  mutate(cpi_vol = slide_dbl(cpi, .f = sd, .before = 30, .after = 30, .complete = T)) %>% 
  drop_na

median_cpi_vol_mo <- 
median_cpi_mo %>% 
  mutate(median_cpi_vol = slide_dbl(median_cpi, .f = sd, .before = 30, .after = 30, .complete = T)) %>% 
  drop_na

median_pce_vol_mo <- 
  median_pce_mo %>% 
  mutate(median_pce_vol = slide_dbl(median_pce, .f = sd, .before = 30, .after = 30, .complete = T)) %>% 
  drop_na

trimmed_cpi_vol_mo <- 
trimmed_cpi_mo %>% 
  mutate(trimmed_cpi_vol = slide_dbl(trimmed_cpi, .f = sd, .before = 30, .after = 30, .complete = T)) %>%
  drop_na
  
df_vol <- 
sbbi_mo %>% 
  right_join(cpi_vol_mo, by = "date") %>% 
  left_join(median_cpi_vol_mo, by = "date") %>% 
  left_join(median_pce_vol_mo, by = "date") %>% 
  left_join(trimmed_cpi_vol_mo, by = "date") %>% 
  mutate(high_median_cpi = ifelse(median_cpi >= 0.0037, 1, 0),
         high_median_pce = ifelse(median_pce >= 0.0037, 1, 0))


# Regressions

# CPI vol

mod_large_cpi_vol <- 
  df_vol %>% 
lm(large ~ cpi_vol , data = .) 

mod_small_cpi_vol <- 
 df_vol %>% 
lm(small ~ cpi_vol, data = .) 

# Median CPI vol

mod_large_median_cpi_vol <- 
  df_vol %>% 
lm(large ~ median_cpi_vol * high_median_cpi, data = .) 

mod_small_median_cpi_vol <- 
 df_vol %>% 
lm(small ~ median_cpi_vol * high_median_cpi, data = .)

# Median PCE vol

mod_large_median_pce_vol <- 
  df_vol %>% 
lm(large ~ median_pce_vol * high_median_pce, data = .) 

mod_small_median_pce_vol <- 
 df_vol %>% 
lm(small ~ median_pce_vol * high_median_pce, data = .)

# Trimmed mean vol

mod_large_trimmed_cpi_vol <- 
  df_vol %>% 
lm(large ~ trimmed_cpi_vol, data = .) 

mod_small_trimmed_cpi_vol <- 
 df_vol %>% 
lm(small ~ trimmed_cpi_vol, data = .)


# Table

rob_se <- 
  list(sqrt(diag(vcovHC(mod_large_cpi_vol))),
      sqrt(diag(vcovHC(mod_large_median_cpi_vol))),
       sqrt(diag(vcovHC(mod_large_median_pce_vol))),
       sqrt(diag(vcovHC(mod_large_trimmed_cpi_vol))),
       sqrt(diag(vcovHC(mod_small_cpi_vol))),
       sqrt(diag(vcovHC(mod_small_median_cpi_vol))),
       sqrt(diag(vcovHC(mod_small_median_pce_vol))),
       sqrt(diag(vcovHC(mod_small_trimmed_cpi_vol))))
       

regressions <- 
  list(mod_large_cpi_vol,
       mod_large_median_cpi_vol,
       mod_large_median_pce_vol,
       mod_large_trimmed_cpi_vol,
       mod_small_cpi_vol,
       mod_small_median_cpi_vol,
       mod_small_median_pce_vol,
       mod_small_trimmed_cpi_vol)

se_list <- 
  list(c("Standard errors", "HC", "HC", "HC", "HC", "HC", "HC", "HC", "HC"))

stargazer(
  regressions,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  digits = 3,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = c("Constant"),
  type = "html",
  out = "regressions_vol.html",
  notes = c("Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted"),
  dep.var.caption = "Dependent variable",
  dep.var.labels = c("Large cap", "Small cap"),
  model.numbers = T,
  add.lines = se_list)

## 
## <table style="text-align:center"><tr><td colspan="9" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="8">Dependent variable</td></tr>
## <tr><td></td><td colspan="8" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="4">Large cap</td><td colspan="4">Small cap</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td><td>(6)</td><td>(7)</td><td>(8)</td></tr>
## <tr><td colspan="9" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">cpi_vol</td><td>-2.146</td><td></td><td></td><td></td><td>-0.564</td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(1.760)</td><td></td><td></td><td></td><td>(2.409)</td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">median_cpi_vol</td><td></td><td>-0.246</td><td></td><td></td><td></td><td>3.102</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(9.882)</td><td></td><td></td><td></td><td>(13.813)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">high_median_cpi</td><td></td><td>0.120<sup>*</sup></td><td></td><td></td><td></td><td>0.120</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.067)</td><td></td><td></td><td></td><td>(0.088)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">median_cpi_vol:high_median_cpi</td><td></td><td>-146.907<sup>*</sup></td><td></td><td></td><td></td><td>-160.313</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(86.620)</td><td></td><td></td><td></td><td>(120.509)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">median_pce_vol</td><td></td><td></td><td>-3.927</td><td></td><td></td><td></td><td>0.791</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(8.716)</td><td></td><td></td><td></td><td>(11.985)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">high_median_pce</td><td></td><td></td><td>-0.045</td><td></td><td></td><td></td><td>-0.021</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(0.079)</td><td></td><td></td><td></td><td>(0.129)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">median_pce_vol:high_median_pce</td><td></td><td></td><td>37.635</td><td></td><td></td><td></td><td>25.919</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(55.056)</td><td></td><td></td><td></td><td>(89.231)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">trimmed_cpi_vol</td><td></td><td></td><td></td><td>1.302</td><td></td><td></td><td></td><td>-0.049</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(8.450)</td><td></td><td></td><td></td><td>(12.016)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="9" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Standard errors</td><td>HC</td><td>HC</td><td>HC</td><td>HC</td><td>HC</td><td>HC</td><td>HC</td><td>HC</td></tr>
## <tr><td style="text-align:left">Observations</td><td>743</td><td>432</td><td>491</td><td>432</td><td>743</td><td>432</td><td>491</td><td>432</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.002</td><td>0.009</td><td>0.004</td><td>0.0001</td><td>0.0001</td><td>0.005</td><td>0.006</td><td>0.00000</td></tr>
## <tr><td colspan="9" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="8" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## <tr><td style="text-align:left"></td><td colspan="8" style="text-align:right">Standard errors are in parentheses are heteroskedasticity adjusted (HC), heteroskedasticity and autocorrelation adjusted (HAC) or IID as noted</td></tr>
## </table>

Figure___: rising median - line chart with shading and table of how SBB have performed

Table __. Median CPI elevated/rising returns

Table (appendix) median CPI elev/rising sd

median_cpi <- 
  c("MEDCPIM158SFRBCLE") %>% 
  tq_get(get = "economic.data", from = "1983-01-01") %>% 
  select(date, median = price) %>%
  mutate(median = ((1 + (median / 100))^(1/12)) - 1,
         median = median + 1,
         median = slide_dbl(median, .f = prod, .before = 11, .complete = T),
         median_cpi = 100*(median - 1)) %>% 
  drop_na %>% 
  select(date, median_cpi)

shading <- 
  median_cpi %>% 
  mutate(rising = ifelse(median_cpi - lag(median_cpi, 6) >= .25 |
                           median_cpi >= 4.5, 1, 0),
         grouping = ifelse(rising == 1 & lead(rising, 1) == 1 |
                             rising == 1 & lead(rising, 2) == 1 |
                             rising == 1 & lead(rising, 3) == 1, 1, 0)) %>% 
  drop_na %>% 
  mutate(year = year(date)) %>% 
  group_by(grouping, year) %>% 
  mutate(begin = first(date), end = last(date)) %>% 
  ungroup %>% 
  filter(grouping == 1) %>% 
  group_by(year) %>% 
  mutate(group_length = length(year)) %>% 
  ungroup

# Identify rising median periods for case_when

shading %>% 
  group_by(begin) %>% 
  filter(row_number() == n()) %>% 
  filter(year < 2023)

## # A tibble: 14 x 8
## # Groups:   begin [14]
##    date       median_cpi rising grouping  year begin      end       
##    <date>          <dbl>  <dbl>    <dbl> <int> <date>     <date>    
##  1 1988-12-01       4.12      1        1  1988 1988-08-01 1988-12-01
##  2 1989-12-01       4.52      1        1  1989 1989-01-01 1989-12-01
##  3 1990-12-01       4.78      1        1  1990 1990-01-01 1990-12-01
##  4 1991-03-01       4.62      1        1  1991 1991-01-01 1991-03-01
##  5 1998-10-01       2.92      1        1  1998 1998-09-01 1998-10-01
##  6 2000-12-01       3.03      1        1  2000 2000-03-01 2000-12-01
##  7 2001-01-01       3.05      1        1  2001 2001-01-01 2001-01-01
##  8 2004-09-01       2.35      1        1  2004 2004-05-01 2004-09-01
##  9 2005-12-01       2.76      1        1  2005 2005-11-01 2005-12-01
## 10 2006-09-01       3.04      1        1  2006 2006-07-01 2006-09-01
## 11 2011-12-01       2.31      1        1  2011 2011-02-01 2011-12-01
## 12 2012-01-01       2.39      1        1  2012 2012-01-01 2012-01-01
## 13 2021-12-01       4.00      1        1  2021 2021-08-01 2021-12-01
## 14 2022-12-01       6.99      1        1  2022 2022-01-01 2022-12-01
## # ... with 1 more variable: group_length <int>

# SBB and median inflation table

read_csv("sbbi_monthly.csv", locale = locale(encoding = "latin1")) %>% 
    select(c(1, 2, 5, 7, 15)) %>% 
    set_colnames(c("date", "large", "small", "lt_gov", "T_bills")) %>% 
  mutate(date = seq(as.Date("1926-01-01"),
                    by = "month",
                    length = length(large))) %>% 
mutate(rising_median =
         case_when(between(date, as.Date("1988-08-01"), as.Date("1991-03-01")) ~ "1988-91",
                   between(date, as.Date("2000-03-01"), as.Date("2001-01-01")) ~ "2000-01",
                   between(date, as.Date("2004-05-01"), as.Date("2004-09-01")) ~ "2004",
                   between(date, as.Date("2011-02-01"), as.Date("2012-01-01")) ~ "2011-12",
                   between(date, as.Date("2021-08-01"), as.Date("2022-12-01")) ~ "2021-22",
                   TRUE ~ "Stable/falling")) %>% 
  filter(year(date) > 1983) %>% 
  mutate(across(large:T_bills, ~ 1 + .)) %>% 
  group_by(rising_median) %>% 
  summarize(months = n(),
            across(large:T_bills, ~ 100*(((last(.)/first(.))^(12/length(large))) - 1))) %>% 
  slice(-n()) %>% 
  add_row(rising_median = "Falling/stable", large = 12.4, small = 17.2, lt_gov = 5.4, T_bills = 3.2) %>% 
  kbl(digits = 1) %>% 
  kable_classic(full_width = F)

rising_median	months	large	small	lt_gov	T_bills
1988-91	32	2.2	3.5	-0.1	-0.1
2000-01	11	-6.2	25.4	-3.8	0.1
2004	5	-0.7	12.4	3.6	0.1
2011-12	12	1.0	1.0	-0.9	0.0
2021-22	17	-6.1	-5.4	-0.8	0.2
Falling/stable	NA	12.4	17.2	5.4	3.2

# Last row of table above (added with add_row) calc

read_csv("sbbi_monthly.csv", locale = locale(encoding = "latin1")) %>% 
    select(c(1, 2, 5, 7, 15)) %>% 
    set_colnames(c("date", "large", "small", "lt_gov", "T_bills")) %>% 
  mutate(date = seq(as.Date("1926-01-01"),
                    by = "month",
                    length = length(large))) %>% 
  filter(year(date) < 2023) %>%
  mutate(rising_median =
         case_when(between(date, as.Date("1988-08-01"), as.Date("1991-03-01")) ~ "1",
                   between(date, as.Date("2000-03-01"), as.Date("2001-01-01")) ~ "1",
                   between(date, as.Date("2004-05-01"), as.Date("2004-09-01")) ~ "1",
                   between(date, as.Date("2011-02-01"), as.Date("2012-01-01")) ~ "1",
                   between(date, as.Date("2021-08-01"), as.Date("2022-12-01")) ~ "1",
                   TRUE ~ "NA")) %>% 
  filter(rising_median != 1) %>% 
  summarize(across(large: T_bills, ~ ((1 + mean(., na.rm = T))^12) - 1))

## # A tibble: 1 x 4
##   large small lt_gov T_bills
##   <dbl> <dbl>  <dbl>   <dbl>
## 1 0.124 0.172 0.0536  0.0317

# SD values for above table

sd_table <- 
  read_csv("sbbi_monthly.csv", locale = locale(encoding = "latin1")) %>% 
    select(c(1, 2, 5, 7, 15)) %>% 
    set_colnames(c("date", "large", "small", "lt_gov", "T_bills")) %>% 
  mutate(date = seq(as.Date("1926-01-01"),
                    by = "month",
                    length = length(large))) %>% 
mutate(rising_median =
         case_when(between(date, as.Date("1988-08-01"), as.Date("1991-03-01")) ~ "1988-91",
                   between(date, as.Date("2000-03-01"), as.Date("2001-01-01")) ~ "2000-01",
                   between(date, as.Date("2004-05-01"), as.Date("2004-09-01")) ~ "2004",
                   between(date, as.Date("2011-02-01"), as.Date("2012-01-01")) ~ "2011-12",
                   between(date, as.Date("2021-08-01"), as.Date("2022-12-01")) ~ "2021-22",
                   TRUE ~ "Stable/falling")) %>% 
  filter(year(date) > 1983) %>% 
  group_by(rising_median) %>% 
  summarize(months = n(),
            across(large:T_bills, ~ sd(.)*sqrt(12)))  %>% 
  add_row(.before = 6, rising_median = "mean", months = NA, large = 0.15, small = 0.22, lt_gov = 0.09, T_bills = .002) %>% 
    kbl(digits = 2) %>% 
  kable_classic(full_width = F)


# Calc means ex. stable/falling
  
# sd_table %>% 
#   select(where(is.numeric)) %>% 
#   slice(-8) %>% 
#   colMeans


## Line chart 

median_cpi %>% 
  ggplot(aes(x = date, y = median_cpi)) +
  geom_rect(aes(xmin = begin, xmax = end, ymin = -Inf, ymax = +Inf), data = shading, fill = "lightgray") +
  geom_rect(aes(xmin = as.Date("1988-01-01"), xmax = as.Date("1991-03-01"), ymin = -Inf, ymax = +Inf), data = shading, fill = "lightgray") +
  # geom_rect(aes(xmin = as.Date("2000-08-01"), xmax = as.Date("2001-06-01"), ymin = -Inf, ymax = +Inf), data = shading, fill = "pink") +
  # geom_rect(aes(xmin = as.Date("2004-10-01"), xmax = as.Date("2005-03-01"), ymin = -Inf, ymax = +Inf), data = shading, fill = "pink") +
  # geom_rect(aes(xmin = as.Date("2011-03-01"), xmax = as.Date("2012-06-01"), ymin = -Inf, ymax = +Inf), data = shading, fill = "pink") +
  geom_rect(aes(xmin = as.Date("2021-10-01"), xmax = as.Date("2023-03-01"), ymin = -Inf, ymax = +Inf), data = shading, fill = "lightgray") +
  geom_line() +
  theme +
  scale_x_date(expand = c(0, 0)) +
  labs(x = NULL, y = "
       Y/y median CPI change  ", caption = "Source: FRED")