Libraries

Note: The user will need to download most of these to run the code in the chunks below.

Note: Theme settings below correspond to those used in the paper’s graphs.

Note: tables created using Stargazer package: Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary StatisticsTables. R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

library(tidyquant) 
library(tidyverse) 
library(knitr) 
library(kableExtra)
library(tsibble)
library(data.table)
library(roll)
library(janitor)
library(fable)
library(mFilter)
library(ggrepel)
library(lmtest)
library(magrittr)
library(sandwich)
library(slider)
library(tseries)
library(stargazer)
theme =  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    plot.title = element_text(
      size = 40,
      lineheight = .9,
      face = "bold",
      hjust = .5,
      vjust = .75,
      family = "serif"
    ),
    plot.subtitle = element_text(hjust = .5, vjust = 1, family = "serif", size = 10),
    plot.caption = element_text(hjust = .5, family = "serif"),
    panel.background = element_rect(fill = "white"),
    axis.line.x = element_line(color = "black"),
    axis.line.y = element_line(color = "black"),
    axis.text = element_text(family = "serif"),
    axis.title = element_text(family = "serif"))

Introduction

Step 1: Create 1960-2022 for “longer” sample period dataframe

Note: Commodity futures data available upon request.

Note the process for detrending PCE involves saving time series objects seprately to apply HP filter. HP filter frequency/lambda value = 129,600 per Hodrick-Prescott, Mankiw, others.

infl =
  c("CPIAUCSL", "CPILFESL", "CPIENGSL", "PCEPI", "PCEPILFE", "UNRATE") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  mutate(across(2:6, ~ 100*((. / lag(., 12)) - 1))) %>% 
  drop_na %>% 
  mutate(month = yearmonth(date),
         cpiaucsl_ma = slide_dbl(cpiaucsl, .f = mean, .before = 18, .after = 18, .complete = T)) %>% 
  select(month, cpiaucsl, cpilfesl, cpiengsl, cpiaucsl_ma, pcepi, pcepilfe, unrate) %>% 
  mutate(unrate_change = unrate - lag(unrate, 12),
         cpi_shock = cpiaucsl - cpilfesl,
         pce_shock = pcepi - pcepilfe)


 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 = yearmonth(date)) %>% 
    select(month, rpce_cyc = cycle)
  
  sp =
    c("^GSPC") %>%
    tq_get(get = "stock.prices", from = "1959-01-01") %>% 
    select(date, adjusted) %>% 
    mutate(month = yearmonth(date)) %>% 
    group_by(month) %>% 
    filter(row_number() == n()) %>% 
    ungroup %>% 
    select(month, adjusted) %>%
    mutate(sp = (adjusted - lag(adjusted, 12)) / lag(adjusted, 12),
           sp_mo = (adjusted - lag(adjusted, 1)) / lag(adjusted, 1)) %>% 
    select(month, sp_mo, sp) %>% 
    drop_na
  
  
  nrou = 
("NROU") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  mutate(month = yearmonth(date)) %>% 
  select(month, nrou = price)
  

  bci =
    read_csv("bci.csv") %>%
    mutate(
      date = mdy(date),
      bcomtr_ret = 100 * ((bcomtr / lag(bcomtr, 12)) - 1),
      bcomtr_ret_mo = 100 * ((bcomtr / lag(bcomtr, 1)) - 1),
      bcom_ret = 100 * ((bcom / lag(bcom, 12)) - 1),
      bcom_ret_mo = 100* ((bcom / lag(bcom, 1)) - 1),
      month = yearmonth(date)) %>%
  drop_na %>%
  select(month, bcomtr_ret, bcomtr_ret_mo, bcom_ret)
  
  
  gsci = 
    read_csv("gsci.csv") %>% 
    select(date, gsci) %>% 
    mutate(date = mdy(date),
           gscitr_ret = 100 * ((gsci / lag(gsci, 12)) - 1),
           month = yearmonth(date)) %>% 
    select(month, gscitr_ret) %>% 
    drop_na

    
  join60 =
  bci %>% 
  left_join(gsci, by = "month") %>% 
  left_join(infl, by = "month") %>% 
  left_join(rpce, by = "month") %>% 
  left_join(sp, by = "month") %>% 
  left_join(nrou, by = "month") %>% 
  fill(nrou, .direction = c("up")) %>% #I fill in missing values of this quarterly series with the next available value
  mutate(tightness = unrate - nrou)

Motivation

Correlation table 1960-2022 (not in paper)

join60 %>% 
  select(where(is.numeric), -c(sp_mo, bcomtr_ret_mo)) %>% 
  cor(use = "complete.obs") %>% 
  stargazer(type = "html", out = "join60 corr.html", digits = 2)
bcomtr_ret bcom_ret gscitr_ret cpiaucsl cpilfesl cpiengsl cpiaucsl_ma pcepi pcepilfe unrate unrate_change cpi_shock pce_shock rpce_cyc sp nrou tightness
bcomtr_ret 1 0.99 0.83 0.51 0.26 0.61 0.44 0.51 0.33 -0.17 -0.17 0.73 0.72 0.21 -0.01 0.35 -0.31
bcom_ret 0.99 1 0.83 0.42 0.15 0.59 0.35 0.42 0.23 -0.20 -0.20 0.75 0.73 0.22 -0.003 0.24 -0.29
gscitr_ret 0.83 0.83 1 0.42 0.17 0.64 0.30 0.41 0.24 -0.24 -0.27 0.70 0.67 0.29 0.07 0.28 -0.34
cpiaucsl 0.51 0.42 0.42 1 0.92 0.70 0.94 0.98 0.91 0.11 0.07 0.46 0.59 0.04 -0.15 0.69 -0.14
cpilfesl 0.26 0.15 0.17 0.92 1 0.43 0.93 0.92 0.96 0.24 0.14 0.09 0.26 -0.05 -0.14 0.76 -0.03
cpiengsl 0.61 0.59 0.64 0.70 0.43 1 0.51 0.66 0.46 -0.11 -0.14 0.83 0.86 0.17 -0.02 0.23 -0.20
cpiaucsl_ma 0.44 0.35 0.30 0.94 0.93 0.51 1 0.95 0.94 0.20 0.13 0.29 0.40 0.01 -0.11 0.78 -0.08
pcepi 0.51 0.42 0.41 0.98 0.92 0.66 0.95 1 0.95 0.17 0.07 0.43 0.54 0.03 -0.14 0.74 -0.10
pcepilfe 0.33 0.23 0.24 0.91 0.96 0.46 0.94 0.95 1 0.29 0.11 0.15 0.26 -0.02 -0.11 0.81 -0.001
unrate -0.17 -0.20 -0.24 0.11 0.24 -0.11 0.20 0.17 0.29 1 0.43 -0.27 -0.26 -0.57 0.02 0.25 0.94
unrate_change -0.17 -0.20 -0.27 0.07 0.14 -0.14 0.13 0.07 0.11 0.43 1 -0.16 -0.10 -0.67 -0.24 0.01 0.44
cpi_shock 0.73 0.75 0.70 0.46 0.09 0.83 0.29 0.43 0.15 -0.27 -0.16 1 0.96 0.20 -0.07 0.03 -0.29
pce_shock 0.72 0.73 0.67 0.59 0.26 0.86 0.40 0.54 0.26 -0.26 -0.10 0.96 1 0.16 -0.14 0.10 -0.30
rpce_cyc 0.21 0.22 0.29 0.04 -0.05 0.17 0.01 0.03 -0.02 -0.57 -0.67 0.20 0.16 1 0.10 0.04 -0.60
sp -0.01 -0.003 0.07 -0.15 -0.14 -0.02 -0.11 -0.14 -0.11 0.02 -0.24 -0.07 -0.14 0.10 1 -0.04 0.03
nrou 0.35 0.24 0.28 0.69 0.76 0.23 0.78 0.74 0.81 0.25 0.01 0.03 0.10 0.04 -0.04 1 -0.11
tightness -0.31 -0.29 -0.34 -0.14 -0.03 -0.20 -0.08 -0.10 -0.001 0.94 0.44 -0.29 -0.30 -0.60 0.03 -0.11 1

Exhibit 1: 1960 to 2022 scatterplots and tables

join60 %>%
  mutate(decade = paste0(10*as.numeric(substr(year(month), 1, 3)), "s")) %>% 
  ggplot(aes(x = cpiaucsl, y = bcomtr_ret)) +
  geom_point(aes(fill = decade), shape = 21, size = 3) +
  geom_smooth(method = "lm", se = F, aes(color = decade)) +
  scale_fill_brewer(palette = "Blues") +
  scale_color_brewer(palette = "Blues") +
  scale_x_continuous(expand = c(0, 0)) +
  theme +
  theme(legend.position = c(1,1), legend.justification = c(1,1), 
        legend.background = element_blank(), legend.title = element_blank(),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        legend.text = element_text(size = 20)) +
  labs(y = "BCOMTR TTM % change", x = "CPI Y/Y % change")

join60 %>%
  mutate(decade = paste0(10*as.numeric(substr(year(month), 1, 3)), "s")) %>% 
  ggplot(aes(x = cpilfesl, y = bcomtr_ret)) +
  geom_point(aes(fill = decade), shape = 21, size = 3) +
  geom_smooth(method = "lm", se = F, aes(color = decade)) +
  scale_fill_brewer(palette = "Blues") +
  scale_color_brewer(palette = "Blues") +
  scale_x_continuous(expand = c(0, 0)) +
  theme +
  theme(legend.position = c(1,1), legend.justification = c(1,1), 
        legend.background = element_blank(), legend.title = element_blank(),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        legend.text = element_text(size = 20)) +
    labs(y = "BCOMTR TTM return %", x = "CPI Y/Y % change")

Regressions for motivation paragraph

# BCI regressions from intro section 

mod1 = 
  join60 %>% 
  lm(bcomtr_ret ~ cpiaucsl, data = .) 

mod2 = 
  join60 %>% 
  lm(bcomtr_ret ~ cpilfesl, data = .) 

mod3 = 
  join60 %>% 
  lm(bcomtr_ret ~ cpiaucsl_ma, data = .)

mod4 = 
  join60 %>% 
  lm(bcomtr_ret ~ cpiaucsl_ma + cpiengsl, data = .) 

mod5 = 
  join60 %>% 
  lm(bcomtr_ret ~ cpilfesl + cpiengsl, data = .) 
   

mod6 = 
  join60 %>% 
  lm(bcomtr_ret ~ pcepi, data = .) 

mod7 = 
  join60 %>% 
  lm(bcomtr_ret ~ pcepilfe, data = .)  

mod8 = 
  join60 %>% 
  lm(bcomtr_ret ~ pcepilfe + cpiengsl, data = .) 


rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))),
              sqrt(diag(vcovHAC(mod6))),
              sqrt(diag(vcovHAC(mod7))),
              sqrt(diag(vcovHAC(mod8))))
stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  mod6,
  mod7,
  mod8,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = c("Constant"),
  type = "html",
  out = "join60 regression table BCOMTR.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted",
            "Source: author's regressions"),
  dep.var.labels = c("Dependent variable: 12-month change in BCOMTR"))
Dependent variable: 12-month change in BCOMTR
(1) (2) (3) (4) (5) (6) (7) (8)
cpiaucsl 4.328***
(1.171)
cpilfesl 2.569*** 0.029
(0.783) (0.909)
cpiaucsl_ma 4.277*** 1.680
(1.479) (1.537)
cpiengsl 1.193*** 1.339*** 1.279***
(0.137) (0.196) (0.162)
pcepi 5.029***
(1.234)
pcepilfe 3.737*** 0.701
(1.089) (1.148)
N 746 746 722 722 746 746 746 746
Adjusted R2 0.257 0.071 0.186 0.382 0.364 0.255 0.107 0.367
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted
Source: author’s regressions 
# Same regressions with BCOM

mod1 =
  join60 %>% 
  lm(bcom_ret ~ cpiaucsl, data = .) 

mod2 =
  join60 %>% 
  lm(bcom_ret ~ cpilfesl, data = .)  

mod3 =
  join60 %>% 
  lm(bcom_ret ~ cpilfesl + cpiengsl, data = .) 

mod4 =
  join60 %>% 
  lm(bcom_ret ~ pcepi, data = .) 

mod5 =
  join60 %>% 
  lm(bcom_ret ~ pcepilfe, data = .)  

mod6 =
  join60 %>% 
  lm(bcom_ret ~ pcepilfe + cpiengsl, data = .) 


rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))),
              sqrt(diag(vcovHAC(mod6))))
stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  mod6,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = c("Constant"),
  type = "html",
  out = "join60 regression table BCOM.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted",
            "Source: author's regressions"),
  dep.var.labels = c("Dependent variable: 12-month change in BCOM"))
Dependent variable: 12-month change in BCOM
(1) (2) (3) (4) (5) (6)
cpiaucsl 3.355***
(1.191)
cpilfesl 1.468* -1.062
(0.794) (0.887)
cpiengsl 1.333*** 1.280***
(0.176) (0.146)
pcepi 3.883***
(1.284)
pcepilfe 2.428** -0.611
(1.156) (1.145)
N 746 746 746 746 746 746
Adjusted R2 0.179 0.026 0.364 0.176 0.052 0.355
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted
Source: author’s regressions 
# Same regressions with GSCI

mod1 =
  join60 %>% 
  lm(gscitr_ret ~ cpiaucsl, data = .) 

mod2 =
  join60 %>% 
  lm(gscitr_ret ~ cpilfesl, data = .)  

mod3 =
  join60 %>% 
  lm(gscitr_ret ~ cpilfesl + cpiengsl, data = .) 

mod4 =
  join60 %>% 
  lm(gscitr_ret ~ pcepi, data = .) 

mod5 =
  join60 %>% 
  lm(gscitr_ret ~ pcepilfe, data = .)  

mod6 =
  join60 %>% 
  lm(gscitr_ret ~ pcepilfe + cpiengsl, data = .) 


rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))),
              sqrt(diag(vcovHAC(mod6))))
stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  mod6,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = c("Constant"),
  type = "html",
  out = "join60 regression table GSCI.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted",
            "Source: author's regressions"),
  dep.var.labels = c("Dependent variable: 12-month change in GSCI"))
Dependent variable: 12-month change in GSCI
(1) (2) (3) (4) (5) (6)
cpiaucsl 3.856**
(1.540)
cpilfesl 1.881* -1.269
(1.087) (1.010)
cpiengsl 1.576*** 1.522***
(0.229) (0.239)
pcepi 4.365***
(1.597)
pcepilfe 2.938** -0.812
(1.444) (1.223)
N 624 624 624 624 624 624
Adjusted R2 0.194 0.034 0.450 0.185 0.062 0.440
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted
Source: author’s regressions

Data

Commodity futures BCOMTR and BCOM returns and SD

# Return and risk

bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         year = as.factor(year(date))) 

bci %>% 
  mutate(across(2:3, ~ (. / lag(., 12)) - 1)) %>% 
  group_by(year) %>% 
  filter(row_number() == n()) %>% 
  ungroup %>% 
  drop_na %>% 
  select(where(is.numeric)) %>% 
  summarize(across(everything(), ~ mean(., na.rm = T))) #change mean to sd for standard deviation
## # A tibble: 1 × 2
##   bcomtr   bcom
##    <dbl>  <dbl>
## 1 0.0988 0.0496
# Correlation with inflation (monthly change)

cpi =
  ("CPIAUCSL") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  select(date, cpi = price) %>% 
  mutate(infl = 100*((cpi/lag(cpi, 1)) - 1))

bci %>% 
  mutate(date = floor_date(date, "month")) %>% 
  left_join(cpi, by = "date") %>% 
  mutate(across(2:3, ~ 100*((. / lag(., 1)) - 1))) %>% 
  filter(year(date) < 2023) %>% 
  drop_na %>% 
  cor.test(~ bcomtr + infl, data = .)
## 
##  Pearson's product-moment correlation
## 
## data:  bcomtr and infl
## t = 4.5465, df = 753, p-value = 6.354e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.09319066 0.23209907
## sample estimates:
##       cor 
## 0.1634549

Exhibit 2: BCI and CPI components other than energy (shown in article)

## med care, used cars, apparel, shelter, trans

cpi =
  c("CPIMEDSL", "CUSR0000SETA02", "CUUR0000SETA01" ,"CPIAPPSL", "CUSR0000SAH1", "CUUR0000SAS4", "CPIENGSL", "CPIUFDSL") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  set_colnames(c("date", "medical", "used_vehicles", "new_vehicles", "apparel", "shelter", "transp", "energy", "food")) %>% 
  mutate(across(2:9, ~100*(. - lag(., 12)) / lag(., 12)),
         month = yearmonth(date)) %>% 
  relocate(month) %>% 
  select(-date)
  
# Exhibit 2: BCI and CPI components

read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         month = yearmonth(date)) %>% 
  drop_na %>% 
  select(month, bcomtr_ret) %>% 
  left_join(cpi, by = "month") %>% 
  select(-month) %>% 
  pivot_longer(-bcomtr_ret) %>% 
  filter(bcomtr_ret < abs(50)) %>% 
  ggplot(aes(x = bcomtr_ret, y = value)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(.~name, ncol = 2, scales = "free_y") +
  theme +
  theme(strip.text = element_text(family = "serif", size = 18, face = "bold"),
        strip.background = element_rect(fill = "lightblue")) +
  labs(x = "BCOMTR TTM return", y = "CPI component inflation %")

Step 2: Create 1983-2022 for “longer” sample period dataframe

Note the process for detrending GDP involves saving time series objects seprately to apply HP filter

Median is weighted by the industries relative importance in the aggregate price index

# PCE for real PCE HP trend

  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(month, rpce_cyc = cycle)

# BCI
  
bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         bcom_ret = 100*(bcom - lag(bcom, 12)) / lag(bcom, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  filter(year(date) < 2023) %>% 
  select(month, bcomtr_ret, bcom_ret) 

# GSCI

 gsci = 
    read_csv("gsci.csv") %>% 
    select(date, gsci) %>% 
    mutate(date = mdy(date),
           gscitr_ret = 100 * ((gsci / lag(gsci, 12)) - 1)) %>% 
    select(month = date, gscitr_ret) %>% 
    drop_na


#SP 500

sp =
  c("^GSPC") %>% 
  tq_get(get = "stock.prices", from = "1980-01-01") %>% 
  select(date, adjusted) %>% 
  mutate(month = floor_date(date, "month"),
         mo = MONTH(date),
         .before = "adjusted") %>% 
  group_by(month) %>% 
  #filter(row_number() == 1 | row_number() == n()) %>% 
  filter(row_number() == n()) %>% 
  ungroup %>% 
  mutate(change_12_mo_sp = 100*(adjusted - lag(adjusted, 12)) / lag(adjusted, 12)) %>% 
  drop_na %>% 
  select(month, change_12_mo_sp)

# FRED macro

infl =
  c("CPIAUCSL", "CPILFESL", "CPIENGSL", "CPIUFDSL", "PCEPI", "PCEPILFE", "AHETPI", "MICH", "EXPINF1YR", "UNRATE") %>%
  tq_get(get = "economic.data", from = "1980-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  select(date, cpiaucsl, cpilfesl, cpiengsl, cpiufdsl, pcepi,  pcepilfe, ahetpi, mich, cleve = expinf1yr, unrate) %>%
  mutate(across(2:8, ~ 100*(. - lag(. ,12)) / lag(., 12)),
         month = floor_date(date, "month"),
         #cpiaucsl_ma = slide_dbl(cpiaucsl, .f = mean, .before = 18, .after = 18, .complete = T)
         ) %>%
  select(-date) %>% 
  relocate(month) %>% 
  drop_na 

# NROU

nrou = 
("NROU") %>% 
  tq_get(get = "economic.data", from = "1980-01-01") %>% 
  mutate(month = floor_date(date, "month")) %>% 
  select(month, nrou = price)
     

# National financial conditions index

nfci =
  c("NFCI") %>%
  tq_get(get = "economic.data", from = "1970-01-01") %>%
  mutate(month = floor_date(date, "month")) %>%
  select(month, price) %>%
  group_by(month) %>%
  filter(row_number() == n()) %>%
  ungroup %>%
  select(month, nfci = price) %>% 
  mutate(nfci_change = nfci - lag(nfci, 12)) %>% 
  select(month, nfci_change)
  

# Median and trimmed 

 mcpi =
   read_csv("mcpi.csv") %>% 
  mutate(date = mdy(date),
         ret = (1 + mcpi),
         mcpi = 100*(roll_prod(ret, width = 12)-1),
         month = floor_date(date, "month")) %>% 
   select(month, mcpi) %>% 
   drop_na
 
 mpce =
   read_csv("mpce.csv") %>% 
     mutate(date = mdy(date),
            ret = (1 + mpce/100),
            mpce = 100*(roll_prod(ret, width = 12)-1),
            month = floor_date(date, "month")) %>% 
     select(month, mpce) %>% 
     drop_na
 
 trimmed =
   read_csv("trimcpi.csv") %>% 
  mutate(date = mdy(date),
         ret = (1 + trim),
         trim = 100*(roll_prod(ret, width = 12)-1),
         month = floor_date(date, "month")) %>% 
   select(month, trim) %>% 
   drop_na
 
    
# Join
 
join83 = 
  bci %>% 
  left_join(gsci, by = "month") %>% 
  left_join(sp, by = "month") %>% 
  left_join(infl, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  left_join(mcpi, by = "month") %>% 
  left_join(mpce, by = "month") %>% 
  left_join(trimmed, by = "month") %>% 
  left_join(rpce, by = "month") %>% 
  left_join(nrou, by = "month") %>% 
  filter(year(month) < 2023) %>%  
  #fill(nrou, .direction = c("up")) %>% 
  rename(sp = change_12_mo_sp) %>% 
  mutate(nrou = case_when(is.na(nrou) & is.na(lead(nrou)) ~ lag(nrou),
                          is.na(nrou) & !is.na(lead(nrou)) ~ lead(nrou),
                          TRUE ~ nrou),
         median_cpi_shock = cpiaucsl - mcpi, 
         median_pce_shock = pcepi - mpce, 
         core_cpi_shock = cpiaucsl - cpilfesl, 
         core_pce_shock = pcepi - pcepilfe, 
         xfe_cpi_shock_median = cpilfesl - mcpi,
         xfe_pce_shock = pcepilfe - mpce,
         xfe_cpi_shock_trim = cpilfesl - trim,
         rpce_cyc = roll_sum(rpce_cyc, width = 12),
         mich_change = mich - lag(mich, 12),
         tightness = nrou - unrate,
         mich_cpi_sur = mcpi - lag(mich, 12),
         mich_pce_sur = mpce - lag(mich, 12),
         cleve_cpi_sur = mcpi - lag(cleve, 12),
         cleve_pce_sur = mpce - lag(cleve, 12))

Exhibit 3: inflation chart (shown in article)

join83 %>%
  select(month, cpiaucsl, cpilfesl, mcpi, trim) %>% 
  drop_na %>% 
  #summarize(across(everything(), sd)) 
  ggplot(aes(x = month)) +
  geom_line(aes(y = mcpi, color = "median"), linewidth = 1.5) +
  geom_line(aes(y = cpilfesl, color = "core"), linewidth = 1.5) +
  geom_line(aes(y = cpiaucsl, color = "headline"), linewidth = 1.5) +
  geom_line(aes(y = trim, color = "trimmed mean"), linewidth = 1.5) +
  scale_x_date(date_labels = "%Y",date_breaks = "3 years", expand = c(0, 0)) +
  scale_color_manual(values = c("steelblue", "azure4", "cornflowerblue", "lightblue"),
                     guide = guide_legend(title = NULL)) +
  theme +
  theme(legend.position = c(0,1), 
        legend.justification = c(0,1), 
        legend.background = element_blank(),
        legend.text = element_text(family = "serif", size = 20),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20)) +
  labs(x = NULL, y = "Inflation % change",
       title = "Trailing-12-month inflation",
       subtitle = "1983 to 2022",
       caption = "Source: St. Louis Federal Reserve")

Exhibit 4: Correlations 1983-2022 chart (shown in article)

join83 %>% 
  select(bcomtr_ret, bcom_ret, gscitr_ret, cpi = cpiaucsl, pcepi, pcepixfe = pcepilfe, cpixfe = cpilfesl, cpiengsl = cpiengsl, median_cpi = mcpi, median_pce = mpce, trim_cpi = trim, tightness, rpce_cyc, nfci_change) %>% 
  cor(use = "complete.obs") %>% 
  stargazer(type = "html", out = "corr83.html", digits = 2)
bcomtr_ret bcom_ret gscitr_ret cpi pcepi pcepixfe cpixfe cpiengsl median_cpi median_pce trim_cpi tightness rpce_cyc nfci_change
bcomtr_ret 1 0.98 0.93 0.66 0.67 0.40 0.27 0.77 0.18 0.23 0.33 0.25 0.15 -0.10
bcom_ret 0.98 1 0.90 0.61 0.61 0.31 0.16 0.79 0.07 0.11 0.24 0.18 0.10 -0.12
gscitr_ret 0.93 0.90 1 0.73 0.72 0.46 0.35 0.79 0.27 0.30 0.42 0.28 0.24 -0.03
cpi 0.66 0.61 0.73 1 0.98 0.80 0.76 0.74 0.70 0.69 0.83 0.34 0.38 0.13
pcepi 0.67 0.61 0.72 0.98 1 0.87 0.77 0.70 0.69 0.72 0.83 0.28 0.36 0.07
pcepixfe 0.40 0.31 0.46 0.80 0.87 1 0.94 0.29 0.83 0.87 0.90 0.24 0.30 -0.002
cpixfe 0.27 0.16 0.35 0.76 0.77 0.94 1 0.15 0.91 0.91 0.94 0.29 0.31 0.07
cpiengsl 0.77 0.79 0.79 0.74 0.70 0.29 0.15 1 0.09 0.10 0.28 0.22 0.26 0.01
median_cpi 0.18 0.07 0.27 0.70 0.69 0.83 0.91 0.09 1 0.95 0.95 0.41 0.39 0.26
median_pce 0.23 0.11 0.30 0.69 0.72 0.87 0.91 0.10 0.95 1 0.92 0.39 0.38 0.19
trim_cpi 0.33 0.24 0.42 0.83 0.83 0.90 0.94 0.28 0.95 0.92 1 0.31 0.36 0.22
tightness 0.25 0.18 0.28 0.34 0.28 0.24 0.29 0.22 0.41 0.39 0.31 1 0.57 0.28
rpce_cyc 0.15 0.10 0.24 0.38 0.36 0.30 0.31 0.26 0.39 0.38 0.36 0.57 1 0.38
nfci_change -0.10 -0.12 -0.03 0.13 0.07 -0.002 0.07 0.01 0.26 0.19 0.22 0.28 0.38 1

Commodity futures and trend inflation

Exhibit 5: regression table and regression results using BCOMTR (shown in article)

#BCI 

mod1 =
  join83 %>% 
  lm(bcomtr_ret ~ cpiaucsl + nfci_change, data = .) 

mod2 =
  join83 %>% 
  lm(bcomtr_ret ~ cpilfesl + nfci_change, data = .) 

mod3 =
  join83 %>% 
  lm(bcomtr_ret ~ mcpi + nfci_change, data = .)

mod4 = 
  join83 %>% 
  lm(bcomtr_ret ~ mcpi + nfci_change + cpiengsl + cpiufdsl, data = .) 

mod5 =
  join83 %>% 
  lm(bcomtr_ret ~ trim + nfci_change, data = .)

mod6 =
  join83 %>% 
  lm(bcomtr_ret ~ trim + nfci_change + cpiengsl + cpiufdsl, data = .) 



rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))),
              sqrt(diag(vcovHAC(mod6))))
  
# Exhibit 5

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  mod6,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "aer",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "join83 regression table.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  covariate.labels = c("CPI", "Core CPI", "Median CPI", "Trimmed-mean CPI", "Fin. conditions change", "Energy CPI", "Food CPI"),
  dep.var.labels = c("Dependent variable: 12-month change in BCOMTR"))
Dependent variable: 12-month change in BCOMTR
(1) (2) (3) (4) (5) (6)
CPI 5.899***
(1.523)
Core CPI 1.731
(2.274)
Median CPI 3.927* 1.612
(2.151) (1.831)
Trimmed-mean CPI 5.786*** 1.714
(1.869) (1.720)
Fin. conditions change -5.030*** -3.860** -5.175 -5.513** -5.980 -5.279**
(1.926) (1.868) (3.712) (2.331) (4.154) (2.338)
Energy CPI 1.113*** 1.089***
(0.109) (0.106)
Food CPI 0.881 0.707
(1.023) (1.063)
Observations 492 492 469 469 469 469
Adjusted R2 0.342 0.039 0.053 0.624 0.137 0.624
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted

Exhibit 5a: regression table and regression results using GSCI

mod1 =
  join83 %>% 
  lm(gscitr_ret ~ cpiaucsl + nfci_change, data = .) 

mod2 =
  join83 %>% 
  lm(gscitr_ret ~ cpilfesl + nfci_change, data = .) 

mod3 =
  join83 %>% 
  lm(gscitr_ret ~ mcpi + nfci_change, data = .)

mod4 = 
  join83 %>% 
  lm(gscitr_ret ~ mcpi + nfci_change + cpiengsl + cpiufdsl, data = .) 

mod5 =
  join83 %>% 
  lm(gscitr_ret ~ trim + nfci_change, data = .)

mod6 =
  join83 %>% 
  lm(gscitr_ret ~ trim + nfci_change + cpiengsl + cpiufdsl, data = .) 



rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))),
              sqrt(diag(vcovHAC(mod6))))
  
# Exhibit 5

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  mod6,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "aer",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "join83 regression table GSCI.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  covariate.labels = c("CPI", "Core CPI", "Median CPI", "Trimmed-mean CPI", "Fin. conditions change", "Energy CPI", "Food CPI"),
  dep.var.labels = c("Dependent variable: 12-month change in GSCI"))
Dependent variable: 12-month change in GSCI
(1) (2) (3) (4) (5) (6)
CPI 9.967***
(1.940)
Core CPI 4.595
(2.955)
Median CPI 7.853*** 4.965***
(2.914) (1.833)
Trimmed-mean CPI 10.391*** 5.304***
(2.283) (1.783)
Fin. conditions change -3.043 -1.000 -5.174 -5.226* -6.238 -4.497
(2.606) (2.575) (4.704) (2.828) (5.153) (2.816)
Energy CPI 1.685*** 1.610***
(0.201) (0.199)
Food CPI 0.819 0.269
(1.416) (1.459)
Observations 492 492 469 469 469 469
Adjusted R2 0.438 0.067 0.083 0.664 0.189 0.666
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted

Exhibit 5b. Regression results using BCOM

mod1 =
  join83 %>% 
  lm(bcom_ret ~ cpiaucsl + nfci_change, data = .) 

mod2 =
  join83 %>% 
  lm(bcom_ret ~ cpilfesl + nfci_change, data = .) 

mod3 = 
  join83 %>% 
  lm(bcom_ret ~ mcpi + nfci_change + cpiengsl + cpiufdsl, data = .) 

mod4 =
  join83 %>% 
  lm(bcom_ret ~ trim + nfci_change + cpiengsl + cpiufdsl, data = .) 

mod5 =
  join83 %>% 
  lm(bcom_ret ~ median_cpi_shock + cpiufdsl + cpiengsl + nfci_change, data = .) 


rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))),
              sqrt(diag(vcovHAC(mod5))))
  

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  mod5,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "join83 regression table BCOM.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 12-month change in BCI"))
Dependent variable: 12-month change in BCI
(1) (2) (3) (4) (5)
cpiaucsl 4.942***
(1.718)
cpilfesl 0.209
(2.480)
mcpi -1.027
(1.744)
trim -0.811
(1.637)
median_cpi_shock 2.210
(2.511)
nfci_change -4.197* -3.265 -5.632** -5.711** -5.046**
(2.388) (2.042) (2.302) (2.268) (2.363)
cpiengsl 1.103*** 1.117*** 0.911***
(0.097) (0.096) (0.252)
cpiufdsl 1.468 1.442 0.884
(0.894) (0.924) (0.568)
N 492 492 469 469 469
Adjusted R2 0.251 0.012 0.653 0.652 0.653
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted

Simple regressions of commodity futures returns for both BCOMTR, BCOM, and GSCI on median PCE (referenced in paper after exhibit 5 discussion)

# Simple regressions of BCOMTR and BCOM on median and trimmed mean for adj r squared


join83 %>% 
  lm(bcomtr_ret ~ mpce, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
mpce 3.458***
(0.450)
Constant -4.428**
(1.724)
Observations 540
R2 0.099
Adjusted R2 0.097
Residual Std. Error 18.433 (df = 538)
F Statistic 59.046*** (df = 1; 538)
Note: p<0.1; p<0.05; p<0.01 
join83 %>% 
  lm(bcom_ret ~ mpce, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcom_ret
mpce 1.619***
(0.440)
Constant -2.781*
(1.685)
Observations 540
R2 0.025
Adjusted R2 0.023
Residual Std. Error 18.017 (df = 538)
F Statistic 13.547*** (df = 1; 538)
Note: p<0.1; p<0.05; p<0.01 
join83 %>% 
  lm(gscitr_ret ~ mpce, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
gscitr_ret
mpce 2.900***
(0.590)
Constant -1.529
(2.260)
Observations 540
R2 0.043
Adjusted R2 0.041
Residual Std. Error 24.168 (df = 538)
F Statistic 24.151*** (df = 1; 538)
Note: p<0.1; p<0.05; p<0.01

Testing Phillips curve components

Commodity futures and shocks

Exhibit 6: Headline XFE shocks (shown in article)

# cpi - median leaves you with food, energy and other shocks. cpi - xfe leaves you with food and enrgy shocks. the difference in these differences leaves you with shocks other than food and energy. headline - median - (headline - core) = core - median

# Exhibit 6: XFE headline shock using median cpi

join83 %>% 
  drop_na %>% 
  mutate(xfe_cpi_shock_ma = slide_dbl(xfe_cpi_shock_median, mean, .before = 12, .after = 0, .complete = T)) %>%
  ggplot(aes(x = month)) +
  geom_line(aes(y = xfe_cpi_shock_ma), color = "azure4", size = 1.5) +
  labs(x = NULL, y = "XFE headline shock %") +
  scale_x_date(date_labels = "%Y", date_breaks = "3 years") +
  scale_y_continuous(breaks = seq(-1, 2, by = 0.5)) +
  theme +
  theme(axis.text.y = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        axis.text.x = element_text(size = 10))

# Shock regressions using BCOMTR 

mod1 =
  join83 %>% 
  lm(bcomtr_ret ~ median_cpi_shock, data = .) #Shocks from everything

mod2 =
  join83 %>% 
  lm(bcomtr_ret ~ core_cpi_shock, data = .) #Shocks from energy/food

mod3 =
  join83 %>% 
  lm(bcomtr_ret ~ xfe_cpi_shock_median, data = .) #Shocks from everything but energy/food using median

mod4 = 
  join83 %>% 
  lm(bcomtr_ret ~ xfe_cpi_shock_trim, data = .) #Shocks from everything but E/F using trimmed mean




rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))


stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = "Constant",
  type = "html",
  out = "join83 CPI shock regression table.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 12-month change in BCI"))
Dependent variable: 12-month change in BCI
(1) (2) (3) (4)
median_cpi_shock 11.468***
(1.031)
core_cpi_shock 11.309***
(1.540)
xfe_cpi_shock_median 9.107**
(4.526)
xfe_cpi_shock_trim -2.723
(5.294)
N 469 492 469 469
R2 0.587 0.469 0.081 0.005
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted 
# Shock regressions using BCOM 

mod1 =
  join83 %>% 
  lm(bcom_ret ~ median_cpi_shock, data = .) #Shocks from everything

mod2 =
  join83 %>% 
  lm(bcom_ret ~ core_cpi_shock, data = .) #Shocks from energy/food

mod3 =
  join83 %>% 
  lm(bcom_ret ~ xfe_cpi_shock_median, data = .) #Shocks from everything but energy/food using median

mod4 = 
  join83 %>% 
  lm(bcom_ret ~ xfe_cpi_shock_trim, data = .) #Shocks from everything but E/F using trimmed mean




rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))


stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = "Constant",
  type = "html",
  out = "join83 CPI shock regression table with BCOM.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 12-month change in BCI"))
Dependent variable: 12-month change in BCI
(1) (2) (3) (4)
median_cpi_shock 11.289***
(0.811)
core_cpi_shock 11.841***
(1.402)
xfe_cpi_shock_median 7.331
(5.369)
xfe_cpi_shock_trim -5.320
(5.790)
N 469 492 469 469
R2 0.608 0.540 0.056 0.019
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted 
# Shock regressions using GSCI


mod1 =
  join83 %>% 
  lm(gscitr_ret ~ median_cpi_shock, data = .) #Shocks from everything

mod2 =
  join83 %>% 
  lm(gscitr_ret ~ core_cpi_shock, data = .) #Shocks from energy/food

mod3 =
  join83 %>% 
  lm(gscitr_ret ~ xfe_cpi_shock_median, data = .) #Shocks from everything but energy/food using median

mod4 = 
  join83 %>% 
  lm(gscitr_ret ~ xfe_cpi_shock_trim, data = .) #Shocks from everything but E/F using trimmed mean




rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))


stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "chi2", "ser", "adj.rsq"),
  omit = "Constant",
  type = "html",
  out = "join83 CPI shock regression table GSCI.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 12-month change in BCI"))
Dependent variable: 12-month change in BCI
(1) (2) (3) (4)
median_cpi_shock 17.192***
(1.890)
core_cpi_shock 16.409***
(2.715)
xfe_cpi_shock_median 14.966**
(5.936)
xfe_cpi_shock_trim -3.068
(7.865)
N 469 492 469 469
R2 0.603 0.470 0.100 0.003
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted

Exhibit 7: shocks from relative price of energy (shown in article)

c("CPIUFDSL", "CPILFESL", "CPIENGSL") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  clean_names %>% 
  mutate(rel_price_food = cpiufdsl/cpilfesl,
         rel_price_energy = cpiengsl/cpilfesl,
         rel_price_food_ma = slide_dbl(rel_price_food, .f = mean, .before = 120, .after = 0, .complete = T),
         rel_price_energy_ma = slide_dbl(rel_price_energy, .f = mean, .before = 120, .after = 0, .complete = T)) %>% 
  filter(year(date) > 1979) %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = rel_price_food_ma, color = "food"), size = 1.5) +
  geom_line(aes(y = rel_price_energy_ma, color = "energy"), size = 1.5) +
  scale_color_manual(values = c("steelblue", "azure4")) +
  scale_x_date(date_labels = "%Y", date_breaks = "3 years", 
               limits = c(as.Date("1980-01-01"), as.Date("2020-01-01")),
               expand = c(0, 0)) +
  scale_y_continuous(breaks = seq(0.5, 1.5, by = 0.05), labels = scales::percent) +
  theme +
  labs(x = NULL, y = "Relative price of commodity") +
  theme(
    legend.position = c(1, 0),
    legend.justification = c(1, 0),
    legend.background = element_blank(),
    legend.title = element_blank(),
    legend.text = element_text(size = 20, family = "serif"),
    axis.text.x = element_text(family = "serif", size = 15),
    axis.text.y = element_text(family = "serif", size = 15),
    axis.ticks.x = element_blank())

Regression for equation 2 (referenced in article) for BCOMTR and GSCI

bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         bcom_ret = 100*(bcom - lag(bcom, 12)) / lag(bcom, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  select(month, bcomtr_ret, bcom_ret) 

# GSCI

 gsci = 
    read_csv("gsci.csv") %>% 
    select(date, gsci) %>% 
    mutate(date = mdy(date),
           gscitr_ret = 100 * ((gsci / lag(gsci, 12)) - 1)) %>% 
    select(month = date, gscitr_ret) %>% 
    drop_na


nfci =
  c("NFCI") %>%
  tq_get(get = "economic.data", from = "1970-01-01") %>%
  mutate(month = floor_date(date, "month")) %>%
  select(month, price) %>%
  group_by(month) %>%
  filter(row_number() == n()) %>%
  ungroup %>%
  select(month, nfci = price) %>% 
  mutate(nfci_change = nfci - lag(nfci, 12)) %>% 
  select(month, nfci_change)

# Regression for equation (2) with BCOMTR

c("CPILFESL", "CPIENGSL") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(rel_price = CPIENGSL/CPILFESL,
         rel_price_ma = slide_dbl(rel_price, .f = mean, .before = 120, .after = 0, .complete = T),
         month = floor_date(date, "month")) %>% 
  select(month, rel_price_ma) %>% 
  left_join(bci, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  drop_na %>% 
  filter(year(month) < 2023) %>% 
  lm(bcomtr_ret ~ rel_price_ma + nfci_change, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
rel_price_ma -78.913***
(8.661)
nfci_change 4.606***
(1.041)
Constant 74.789***
(7.040)
Observations 612
R2 0.151
Adjusted R2 0.148
Residual Std. Error 24.323 (df = 609)
F Statistic 53.993*** (df = 2; 609)
Note: p<0.1; p<0.05; p<0.01 
  #coeftest(vcov. = vcovHAC)

# Regression for equation (2) with GSCITR

c("CPILFESL", "CPIENGSL") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(rel_price = CPIENGSL/CPILFESL,
         rel_price_ma = slide_dbl(rel_price, .f = mean, .before = 120, .after = 0, .complete = T),
         month = floor_date(date, "month")) %>% 
  select(month, rel_price_ma) %>% 
  left_join(gsci, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  drop_na %>% 
  filter(year(month) < 2023) %>% 
  lm(gscitr_ret ~ rel_price_ma + nfci_change, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
gscitr_ret
rel_price_ma -52.721***
(9.021)
nfci_change 3.090***
(1.085)
Constant 52.625***
(7.333)
Observations 612
R2 0.068
Adjusted R2 0.065
Residual Std. Error 25.336 (df = 609)
F Statistic 22.247*** (df = 2; 609)
Note: p<0.1; p<0.05; p<0.01 
  #coeftest(vcov. = vcovHAC)


# Using PCEPI instead of CPI

c("DNRGRC1M027SBEA", "PCEPILFE") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(rel_price = DNRGRC1M027SBEA/PCEPILFE,
         rel_price_ma = slide_dbl(rel_price, .f = mean, .before = 120, .after = 0, .complete = T),
         month = floor_date(date, "month")) %>% 
  select(month, rel_price_ma) %>% 
  left_join(bci, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  drop_na %>% 
  lm(bcomtr_ret ~ rel_price_ma + nfci_change, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
rel_price_ma -9.023***
(0.776)
nfci_change 5.336***
(1.001)
Constant 45.905***
(3.128)
Observations 614
R2 0.210
Adjusted R2 0.207
Residual Std. Error 23.432 (df = 611)
F Statistic 81.020*** (df = 2; 611)
Note: p<0.1; p<0.05; p<0.01 
  #coeftest(vcov. = vcovHAC)

Commodity futures and shocks: relative price of food

# Regression

bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  select(month, bcomtr_ret) 

nfci =
  c("NFCI") %>%
  tq_get(get = "economic.data", from = "1970-01-01") %>%
  mutate(month = floor_date(date, "month")) %>%
  select(month, price) %>%
  group_by(month) %>%
  filter(row_number() == n()) %>%
  ungroup %>%
  select(month, nfci = price) %>% 
  mutate(nfci_change = nfci - lag(nfci, 12)) %>% 
  select(month, nfci_change)

c("CPILFESL", "CPIUFDSL") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(rel_price = CPIUFDSL/CPILFESL,
         rel_price_ma = slide_dbl(rel_price, .f = mean, .before = 120, .after = 0, .complete = T),
         month = floor_date(date, "month")) %>% 
  select(month, rel_price_ma) %>% 
  left_join(bci, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  drop_na %>% 
  lm(bcomtr_ret ~ rel_price_ma + nfci_change, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
rel_price_ma 14.036
(24.191)
nfci_change 5.214***
(1.108)
Constant -2.646
(23.879)
Observations 615
R2 0.035
Adjusted R2 0.032
Residual Std. Error 25.886 (df = 612)
F Statistic 11.119*** (df = 2; 612)
Note: p<0.1; p<0.05; p<0.01 
  #coeftest(vcov. = vcovHAC)


# Using PCEPI instead of CPI

c("DNRGRC1M027SBEA", "PCEPILFE") %>% 
  tq_get(get = "economic.data", from = "1960-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(rel_price = DNRGRC1M027SBEA/PCEPILFE,
         rel_price_ma = slide_dbl(rel_price, .f = mean, .before = 120, .after = 0, .complete = T),
         month = floor_date(date, "month")) %>% 
  select(month, rel_price_ma) %>% 
  left_join(bci, by = "month") %>% 
  left_join(nfci, by = "month") %>% 
  drop_na %>% 
  lm(bcomtr_ret ~ rel_price_ma + nfci_change, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
rel_price_ma -9.023***
(0.776)
nfci_change 5.336***
(1.001)
Constant 45.905***
(3.128)
Observations 614
R2 0.210
Adjusted R2 0.207
Residual Std. Error 23.432 (df = 611)
F Statistic 81.020*** (df = 2; 611)
Note: p<0.1; p<0.05; p<0.01 
  #coeftest(vcov. = vcovHAC)

Commodity futures and unexpected inflation

bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         bcom_ret = 100*(bcom - lag(bcom, 12)) / lag(bcom, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  select(month, bcomtr_ret, bcom_ret) 

# GSCI

 gsci = 
    read_csv("gsci.csv") %>% 
    select(date, gsci) %>% 
    mutate(date = mdy(date),
           gscitr_ret = 100 * ((gsci / lag(gsci, 12)) - 1)) %>% 
    select(month = date, gscitr_ret) %>% 
    drop_na

joinexp1yr = 
c("CPIAUCSL", "CPIENGSL", "CPIUFDSL", "CPILFESL", "PCEPI", "PCEPILFE", "MICH", "EXPINF1YR") %>% 
  tq_get(get = "economic.data", from = "1978-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(across(2:7, ~ 100*(. - lag(., 12)) / lag(., 12)), 
         month = floor_date(date, "month")) %>% 
  select(-date) %>% 
  left_join(bci, by = "month") %>% 
  left_join(gsci, by = "month") %>% 
  relocate(month) %>% 
  mutate(mich_cpi_sur = CPIAUCSL - lag(MICH, 12),
         mich_pce_sur = PCEPI - lag(MICH, 12),
         cleve_cpi_sur = CPIAUCSL - lag(EXPINF1YR, 12),
         cleve_pce_sur = PCEPI - lag(EXPINF1YR, 12), 
         rel_energy = CPIENGSL/CPILFESL,
         rel_energy_ma = slide_dbl(rel_energy, .f = mean, .before = 12, .after = 0, .complete = T))


# Expectations regressions using BCOMTR

mod1 =
  joinexp1yr %>% 
  lm(bcomtr_ret ~ mich_cpi_sur , data = .) 

mod2 = 
  joinexp1yr %>% 
  lm(bcomtr_ret ~ mich_cpi_sur + CPIENGSL, data = .) 


mod3 = 
  joinexp1yr %>% 
  lm(bcomtr_ret ~ cleve_cpi_sur, data = .) 

mod4 = 
  joinexp1yr %>% 
  lm(bcomtr_ret ~ cleve_cpi_sur + CPIENGSL, data = .) 



rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "Expectations 1y regression table BCOMTR.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 60-month compound annual change in BCOMTR"))
Dependent variable: 60-month compound annual change in BCOMTR
(1) (2) (3) (4)
mich_cpi_sur 7.789*** 3.280***
(0.755) (1.060)
CPIENGSL 0.803*** 1.403***
(0.169) (0.185)
cleve_cpi_sur 5.572*** -1.936*
(1.195) (1.122)
N 531 531 483 483
Adjusted R2 0.529 0.607 0.309 0.583
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted 
# Expectations regressions using BCOM 


mod1 =
  joinexp1yr %>% 
  lm(bcom_ret ~ mich_cpi_sur + CPIENGSL + CPIUFDSL , data = .) 

mod2 = 
  joinexp1yr %>% 
  lm(bcom_ret ~ mich_cpi_sur + CPIENGSL, data = .) 


mod3 = 
  joinexp1yr %>% 
  lm(bcom_ret ~ cleve_cpi_sur, data = .) 

mod4 = 
  joinexp1yr %>% 
  lm(bcom_ret ~ cleve_cpi_sur + CPIENGSL, data = .) 



rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "Expectations 1y regression table using BCOM.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 60-month compound annual change in BCOM"))
Dependent variable: 60-month compound annual change in BCOM
(1) (2) (3) (4)
mich_cpi_sur 3.234** 1.541
(1.349) (0.952)
CPIENGSL 0.836*** 0.928*** 1.236***
(0.183) (0.164) (0.162)
CPIUFDSL -1.337*
(0.714)
cleve_cpi_sur 6.001*** -0.613
(1.034) (1.025)
N 531 531 483 483
Adjusted R2 0.583 0.566 0.387 0.616
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted 
# Expectations regressions using GSCI

mod1 =
  joinexp1yr %>% 
  lm(gscitr_ret ~ mich_cpi_sur , data = .) 

mod2 = 
  joinexp1yr %>% 
  lm(gscitr_ret ~ mich_cpi_sur + CPIENGSL, data = .) 


mod3 = 
  joinexp1yr %>% 
  lm(gscitr_ret ~ cleve_cpi_sur, data = .) 

mod4 = 
  joinexp1yr %>% 
  lm(gscitr_ret ~ cleve_cpi_sur + CPIENGSL, data = .) 



rob_se = list(sqrt(diag(vcovHAC(mod1))),
              sqrt(diag(vcovHAC(mod2))),
              sqrt(diag(vcovHAC(mod3))),
              sqrt(diag(vcovHAC(mod4))))

stargazer(
  mod1,
  mod2,
  mod3,
  mod4,
  se = rob_se,
  ci = F,
  ci.level = 0.95,
  style = "qje",
  digits = 3,
  intercept.bottom = FALSE,
  omit.stat = c("f", "rsq", "chi2", "ser"),
  omit = "Constant",
  type = "html",
  out = "Expectations 1y regression table GSCITR.html",
  notes = c("Standard errors are heteroscedasticity and autocorrelation adjusted"),
  dep.var.labels = c("Dependent variable: 60-month compound annual change in GSCITR"))
Dependent variable: 60-month compound annual change in GSCITR
(1) (2) (3) (4)
mich_cpi_sur 9.350*** 4.006**
(1.684) (1.987)
CPIENGSL 0.952*** 1.922***
(0.363) (0.271)
cleve_cpi_sur 8.980*** -1.304
(1.683) (1.186)
N 531 531 483 483
Adjusted R2 0.473 0.541 0.375 0.615
Notes: ***Significant at the 1 percent level.
**Significant at the 5 percent level.
*Significant at the 10 percent level.
Standard errors are heteroscedasticity and autocorrelation adjusted

Exhibit 8: Expectations/surprises 5 year and regression results (shown in article)

Commodity futures and economic slack

HP filter lambda = 1,600 (quarterly data) per Hodrick and Prescott, Mankiw, others

# GDP HP filter

  rgdpts =
  c("GDPC1") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>%
  select(price) %>% 
  ts(start = 1959, frequency = 4)
  
  filter =
  rgdpts %>% 
  hpfilter(drift = F, type = c("lambda"), freq = 1600) 
  
  rgdp =
    tibble(trend = as.vector(filter$trend),
       cycle = as.vector(filter$cycle),
       date = seq(as.Date("1959-01-01"),
                              by = "quarter", 
                              length = length(filter$trend))) %>% 
    mutate(qtr = floor_date(date, "quarter")) %>% 
    select(qtr, trend, cycle)
  
  
  # CPI energy
  
  cpiengsl = 
    c("CPIENGSL") %>% 
    tq_get(get = "economic.data", from = "1959-01-01") %>% 
    mutate(qtr = floor_date(date, "quarter")) %>% 
    group_by(qtr) %>% 
    filter(row_number() == n()) %>% 
    ungroup %>% 
    mutate(cpiengsl_qtr = 100*(price - lag(price, 4)) / lag(price, 4)) %>% 
    select(qtr, cpiengsl_qtr)
  
  # Regression using BCOMTR
  
  read_csv("bci.csv") %>% 
    mutate(date = mdy(date),
           qtr = floor_date(date, "quarter")) %>% 
    group_by(qtr) %>% 
    filter(row_number() == n()) %>% 
    select(date, qtr, bcomtr) %>% 
    ungroup %>% 
    mutate(bcomtr_ret_qtr = 100*(bcomtr - lag(bcomtr, 4)) / lag(bcomtr, 4)) %>% 
    drop_na %>% 
    left_join(rgdp, by = "qtr") %>% 
    left_join(cpiengsl, by = "qtr") %>% 
    filter(year(date) < 2023) %>% 
    #lm(bcomtr_ret_qtr ~ cycle + cpiengsl_qtr, data = .) %>% 
    lm(bcomtr_ret_qtr ~ cycle, data = .) %>% 
    #stargazer(type = "html")
    coeftest(vcov. = vcovHAC)

t test of coefficients:

         Estimate Std. Error t value  Pr(>|t|)

(Intercept) 10.373025 3.112878 3.3323 0.0009939 ** cycle 0.032648 0.014728 2.2167 0.0275557
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

   # Regression using GSCI
  
  read_csv("gsci.csv") %>% 
    mutate(date = mdy(date),
           qtr = floor_date(date, "quarter")) %>% 
    group_by(qtr) %>% 
    filter(row_number() == n()) %>% 
    select(date, qtr, gsci) %>% 
    ungroup %>% 
    mutate(gscitr_ret_qtr = 100*(gsci - lag(gsci, 4)) / lag(gsci, 4)) %>% 
    drop_na %>% 
    left_join(rgdp, by = "qtr") %>% 
    left_join(cpiengsl, by = "qtr") %>% 
    filter(year(date) < 2023) %>% 
    #lm(gscitr_ret_qtr ~ cycle + cpiengsl_qtr, data = .) %>% 
    lm(gscitr_ret_qtr ~ cycle, data = .) %>% 
    # stargazer(type = "html")
    coeftest(vcov. = vcovHAC)

t test of coefficients:

         Estimate Std. Error t value Pr(>|t|)

(Intercept) 10.074873 2.996871 3.3618 0.000923 * cycle 0.046121 0.014800 3.1163 0.002093 — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Economic overheating: NROU - regression results (referenced in article)

unemp =
  c("NROU", "UNRATE") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  pivot_wider(names_from = symbol, values_from = price) %>% 
  mutate(excess_unemp = NROU - UNRATE,
         qtr = floor_date(date, "quarter")) %>% 
  select(qtr, excess_unemp)


# Regression


# CPI energy
  
  cpiengsl = 
    c("CPIENGSL") %>% 
    tq_get(get = "economic.data", from = "1959-01-01") %>% 
    mutate(qtr = floor_date(date, "quarter")) %>% 
    group_by(qtr) %>% 
    filter(row_number() == n()) %>% 
    ungroup %>% 
    mutate(cpiengsl_qtr = 100*(price - lag(price, 4)) / lag(price, 4)) %>% 
    select(qtr, cpiengsl_qtr)
  
  
  # Regression using BCOMTR
  
  read_csv("bci.csv") %>% 
    mutate(date = mdy(date),
           qtr = floor_date(date, "quarter")) %>% 
    group_by(qtr) %>% 
    filter(row_number() == n()) %>% 
    select(qtr, bcom, bcomtr) %>% 
    ungroup %>% 
    mutate(bcomtr_ret_qtr = 100*(bcomtr - lag(bcomtr, 4)) / lag(bcomtr, 4)) %>% 
    drop_na %>% 
    left_join(unemp, by = "qtr") %>% 
    left_join(cpiengsl, by = "qtr") %>% 
    drop_na %>% #to get rid of excess quarterly values in unemp dataset 
    lm(bcomtr_ret_qtr ~ excess_unemp + cpiengsl_qtr, data = .) %>% 
    stargazer(type = "html")
Dependent variable:
bcomtr_ret_qtr
excess_unemp 2.019***
(0.685)
cpiengsl_qtr 1.273***
(0.107)
Constant 5.198***
(1.349)
Observations 249
R2 0.399
Adjusted R2 0.394
Residual Std. Error 18.480 (df = 246)
F Statistic 81.530*** (df = 2; 246)
Note: p<0.1; p<0.05; p<0.01 
    #coeftest(vcov. = vcovHAC)

History paints a rosy picture

Data for exhibits 9 and 10: When trend has outpaced headline (shown in article)

  mpce =
  read_csv("mpce.csv") %>%  
  mutate(mpce_12_mo = 100*((roll_prod(x = (1 + mpce/100), width = 12))-1),
         date = mdy(date),
         month = floor_date(date, "month")) %>% 
  select(month, mpce_12_mo) 

pce = 
  c("PCEPI") %>% 
  tq_get(get = "economic.data", from = "1977-01-01") %>% 
  mutate(pce_12_mo = 100*((price - lag(price, 12)) / lag(price ,12)),
         month  = floor_date(date, "month")) %>% 
  select(month, pce_12_mo) 

sp =
    c("^GSPC") %>%
    tq_get(get = "stock.prices", from = "1959-01-01") %>% 
    select(date, adjusted) %>% 
    mutate(month = floor_date(date, "month")) %>% 
    group_by(month) %>% 
    filter(row_number() == n()) %>% 
    ungroup %>% 
    select(month, adjusted) %>%
    mutate(sp = 100*(adjusted - lag(adjusted, 12)) / lag(adjusted, 12),
           sp_mo = (adjusted - lag(adjusted, 1)) / lag(adjusted, 1)) %>% 
    select(month, sp) %>% 
    drop_na

 bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         bcom_ret = 100*(bcom - lag(bcom, 12)) / lag(bcom, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  select(month, bcomtr_ret, bcom_ret)

Exhibit 9. Investment performance during rising trend inflation episodes

 mpce %>%
   left_join(pce, by = "month") %>%
   left_join(sp, by = "month") %>%
   left_join(bci, by = "month") %>%
   drop_na %>%
   mutate(mpce_less_pce = mpce_12_mo - pce_12_mo,
          bci_less_sp = bcomtr_ret - sp,
          dummy = case_when(mpce_12_mo > lag(mpce_12_mo, 12) & pce_12_mo < lag(pce_12_mo, 12) ~ 1,
                            TRUE ~ 0),
          year = year(month)) %>% 
   filter(dummy == 1) %>% 
   group_by(year) %>% 
   summarize(n = n(),
             mean_sp = mean(sp),
             mean_bcomtr = mean(bcomtr_ret),
             mean_bcom = mean(bcom_ret)) %>% 
   pivot_longer(-c(year, n)) %>% 
   ungroup %>% 
   filter(n >= 3) %>% 
   ggplot(aes(x = factor(year), y = value, fill = name)) +
   geom_col(aes(fill = name), position = "dodge") +
   geom_text(aes(y = value, label = paste0(round(value, 1), "%")), position = position_dodge(width = .9), vjust = -.5, family = "serif", size = 4, font = "bold") +
   theme +
   scale_fill_manual(values = c("azure4", "azure3", "steelblue")) +
   labs(x = NULL, y = "Return") +
   theme(
     legend.position = c(0, 0),
     legend.justification = c(0, 0),
     legend.background = element_blank(),
     legend.title = element_blank(),
     legend.text = element_text(size = 20, family = "serif"),
     axis.text.x = element_text(size = 15),
     axis.text.y = element_text(size = 15)
   )

Exhibit 10. BCOMTR excess return (BCOMTR – S&P 500), vertical axis, versus difference between median and headline PCE, horizontal axis

 mpce %>%
   left_join(pce, by = "month") %>%
   left_join(sp, by = "month") %>%
   left_join(bci, by = "month") %>%
   drop_na %>%
   mutate(mpce_less_pce = mpce_12_mo - pce_12_mo,
          bci_less_sp = bcomtr_ret - sp) %>% 
   ggplot(aes(x = mpce_less_pce, y = bci_less_sp)) +
   geom_hline(yintercept = 0, linetype = 2) +
   geom_vline(xintercept = 0, linetype = 2) +
   geom_point(shape = 21, aes(fill = pce_12_mo), size = 4) +
   geom_text_repel(aes(label = year(month)), vjust = -1, max.overlaps = 1, size = 5) +
   scale_fill_gradient2(low = "cornflowerblue", high = "azure4", mid = "white", midpoint = 5, name = "PCE 12-mo\n% change") +
   theme +
   theme(
     legend.position = c(1, 1),
     legend.justification = c(1, 1),
     legend.background = element_blank(),
     legend.title = element_text(family = "serif"),
     legend.text = element_text(size = 15, family = "serif"),
     axis.text.x = element_text(size = 15),
     axis.text.y = element_text(size = 15)) +
   labs(x = "median PCE less headline PCE", y = "BCI less S&P 500")

Misc

When 37 month centered MA CPI is rising (referenced in article)

Centered 37 mo CPI a decent measure of trend per Cecchetti, Ball, others

cpi =
  c("CPIAUCSL") %>% 
  tq_get(get = "economic.data", from = "1959-01-01") %>% 
  mutate(cpi = 100*(price - lag(price, 12)) / lag(price, 12),
         cpi_ma = slide_dbl(cpi, .f = mean, .before = 30, .after = 30),
         month = floor_date(date, "month")) %>% 
  select(month, cpi, cpi_ma) %>% 
  drop_na

bci =
  read_csv("bci.csv") %>% 
  mutate(date = mdy(date),
         bcomtr_ret = 100*(bcomtr - lag(bcomtr, 12)) / lag(bcomtr, 12),
         bcom_ret = 100*(bcom - lag(bcom, 12)) / lag(bcom, 12),
         month = floor_date(date, "month")) %>% 
  drop_na %>% 
  select(month, bcomtr_ret, bcom_ret) 


cpi %>% 
  left_join(bci, by = "month") %>% 
  lm(bcomtr_ret ~ cpi, data = .) %>% 
  stargazer(type = "html")
Dependent variable:
bcomtr_ret
cpi 4.356***
(0.275)
Constant -6.318***
(1.322)
Observations 729
R2 0.257
Adjusted R2 0.256
Residual Std. Error 21.136 (df = 727)
F Statistic 251.490*** (df = 1; 727)
Note: p<0.1; p<0.05; p<0.01

Footnote calculation: Which core/trend measure has lowest sd of difference with MA CPI?

join83 %>%
  mutate(cpiaucsl_ma = slide_dbl(cpiaucsl, .f = mean, .before = 18, .after = 18, .complete = T)) %>%
  drop_na %>%
  select(month, mcpi, trim, cpilfesl, cpiaucsl_ma) %>%
  mutate(across(2:4, ~. - cpiaucsl_ma)) %>% #The difference in each period between measure and the true trend in inflation 
  summarize(across(2:4, sd)) #The SD of the difference between each trend measure and the true trend. trimmed mean wins
## # A tibble: 1 × 3
##    mcpi  trim cpilfesl
##   <dbl> <dbl>    <dbl>
## 1 0.627 0.518    0.680