Load libraries and personal ggplot theme

library(tidyquant)
library(tidyverse)
library(broom)
library(knitr)
library(kableExtra)
library(tsibble)
library(ggcorrplot)
library(ggfittext)
library(psych)
library(data.table)
library(ivreg)
library(janitor)
library(readxl)
library(fable)
library(mFilter)
library(ggrepel)
library(lmtest)
library(magrittr)
library(sandwich)
library(car)
library(slider)
library(tsibble)
library(feasts)
library(tseries)
library(ggfittext) # for geom_fit_text
library(ggalt) # for geom_encircle
library(stargazer) # for markdown printing: {r results = "asis", message=FALSE, warning=FALSE}
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"))

Load data

sbbi <- 
  read_csv("sbbi_july.csv", locale = locale(encoding = "latin1"), col_names = F) %>% 
        select(c(1, 2, 5)) %>% 
        set_colnames(c("date", "large", "small")) %>% 
        mutate(date = seq(as.Date("1926-01-01"),
                          by = "month",
                          length = length(large))) %>% 
  drop_na


df <- 
  read_csv("long_rate.csv") %>% 
  mutate(date = seq(as.Date("1871-01-01"),
                    by = "month",
                    length = length(rate))) %>% 
  left_join(sbbi, by = "date") %>% 
  drop_na %>% 
  mutate(rate_chg = rate - lag(rate, 1)) %>% 
  select(date, rate, rate_chg, small, large) %>% 
  mutate(index = row_number())

Correlations (stocks ~ rates and lag(rates))

df %>% 
  #filter(rate_chg > 0) %>% 
  cor.test(~ large + rate_chg, data = .)

## 
##  Pearson's product-moment correlation
## 
## data:  large and rate_chg
## t = -3.3805, df = 1177, p-value = 0.0007472
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.15428867 -0.04119999
## sample estimates:
##         cor 
## -0.09806089

df %>% 
  #filter(rate_chg > 0) %>% 
  cor.test(~ large + rate_chg, data = .)

## 
##  Pearson's product-moment correlation
## 
## data:  large and rate_chg
## t = -3.3805, df = 1177, p-value = 0.0007472
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.15428867 -0.04119999
## sample estimates:
##         cor 
## -0.09806089

df %>% 
  #filter(rate_chg > 0) %>% 
  cor.test(~ small + lag(rate_chg), data = .)

## 
##  Pearson's product-moment correlation
## 
## data:  small and lag(rate_chg)
## t = -1.3937, df = 1176, p-value = 0.1637
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.09749824  0.01654559
## sample estimates:
##         cor 
## -0.04060858

df %>% 
  #filter(rate_chg > 0) %>% 
  cor.test(~ large + lag(rate_chg), data = .)

## 
##  Pearson's product-moment correlation
## 
## data:  large and lag(rate_chg)
## t = -2.4699, df = 1176, p-value = 0.01366
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.12842618 -0.01478218
## sample estimates:
##         cor 
## -0.07183732

Chart 1 left and right

df %>% 
  ggplot(aes(x = rate_chg, y = small)) +
  geom_point(shape = 21, fill = "dodgerblue", size = 2) +
  theme +
  labs(x = "Long-bond rate change", y = "Small-stock index return", 
       caption = "Source: CFA Institute, Robert Shiller Data")

df %>% 
  ggplot(aes(x = rate_chg, y = large)) +
  geom_point(shape = 21, fill = "dodgerblue", size = 2) +
  theme +
  labs(x = "Long-bond rate change", y = "Large-stock index return", 
       caption = "Source: CFA Institute, Robert Shiller Data")

Table 1 and Table 2(comment in filter for Table 2): by decade

df %>% 
  #filter(rate_chg > 0) %>% 
  mutate(decade = paste0(10*as.numeric(substr(year(date), 1, 3)), "s")) %>% 
  group_by(decade) %>% 
  reframe(small_rate_cor = cor(rate_chg, small, use = c("complete.obs")),
          large_rate_cor = cor(rate_chg, large, use = c("complete.obs"))) %>% 
  kbl(digits = 3) %>% 
  kable_classic(full_width = F,
              font_size = 20,
              html_font = "serif") %>% 
              footnote(general = "\nSource: Robert Shiller data, CFA Institute",
                       footnote_as_chunk = T,
                       general_title = "")

decade	small_rate_cor	large_rate_cor
1920s	0.325	0.271
1930s	-0.133	-0.151
1940s	-0.047	-0.071
1950s	0.063	0.089
1960s	-0.207	-0.254
1970s	-0.272	-0.398
1980s	-0.228	-0.262
1990s	-0.111	-0.263
2000s	0.088	0.103
2010s	0.377	0.306
2020s	-0.009	-0.190
Source: Robert Shiller data, CFA Institute

Chart 2: Rolling correlations

# Small

window <- 60

results_vector <- 
  tibble(values = rep(0, length(df$index) - window))

for(i in 1:length(results_vector$values)){
  
  cor <- 
    df %>% 
    filter(between(index, first(index) + i - 1, first(index) + window + i - 1)) %>% 
    cor.test(~ small + rate_chg, data = .)
  
  results_vector[i,1] <- cor[4] 
}


results_vector %>% 
  mutate(date = df$date[df$index > window]) %>% 
  ggplot(aes(x = date, y = values)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = 1) +
  theme +
  scale_x_date(date_breaks = "5 years", date_labels = "%Y") +
  labs(title = paste0(window, "-month correlation: Small stocks and long-rate change"),
       x = NULL, 
       y = "Correlation",
       caption = "Source: CFAI Institute and Robert Shiller data") +
  theme(axis.text.x = element_text(size = 8))

# Large

window <- 60

results_vector <- 
  tibble(values = rep(0, length(df$index) - window))

for(i in 1:length(results_vector$values)){
  
  cor <- 
    df %>% 
    filter(between(index, first(index) + i - 1, first(index) + window + i - 1)) %>% 
    cor.test(~ large + rate_chg, data = .)
  
  results_vector[i,1] <- cor[4] 
}


results_vector %>% 
  mutate(date = df$date[df$index > window]) %>% 
  ggplot(aes(x = date, y = values)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = 1) + 
  theme +
  scale_x_date(date_breaks = "5 years", date_labels = "%Y") +
  labs(title = paste0(window, "-month correlation: Large stocks and long-rate change"),
       x = NULL, 
       y = "Correlation",
       caption = "Source: CFAI Institute and Robert Shiller data") +
  theme(axis.text.x = element_text(size = 8))

Chart 3: Partial correlations

# Extract residuals

mod <- 
  df %>% 
  lm(small ~ large, data = .) %>% 
  augment %>% 
  select(.resid) %>% 
  mutate(date = df$date,
         rate_chg = df$rate_chg, 
         index = row_number()) 

# Partial corr coeff

mod %>% 
  cor.test(~ .resid + rate_chg, data = .)

## 
##  Pearson's product-moment correlation
## 
## data:  .resid and rate_chg
## t = 1.9467, df = 1177, p-value = 0.05181
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.000442515  0.113375350
## sample estimates:
##        cor 
## 0.05665048

# Rolling partial corr coeff

window <- 60

results_vector <- 
  tibble(values = rep(0, length(df$index) - window))

for(i in 1:length(results_vector$values)){
  
  cor <- 
    mod %>% 
    filter(between(index, first(index) + i - 1, first(index) + window + i - 1)) %>% 
    cor.test(~ .resid + rate_chg, data = .)
  
  results_vector[i,1] <- cor[4] 
}


results_vector %>% 
  mutate(date = mod$date[mod$index > window]) %>% 
  ggplot(aes(x = date, y = values)) +
  geom_line() +
  theme +
  scale_x_date(date_breaks = "5 years", date_labels = "%Y") +
  labs(title = paste0(window, "-month partial correlation: Small stocks and long-rate change"),
       x = NULL, 
       y = "Correlation",
       caption = "Source: CFAI Institute and Robert Shiller data") +
  theme(axis.text.x = element_text(size = 8))

Table 3: Fed tightenings

df %>% 
  mutate(
  tightening = 
    case_when(between(date, as.Date("1965-09-01"), as.Date("1966-11-01")) ~ "1965-6",
              between(date, as.Date("1967-12-01"), as.Date("1969-08-01")) ~ "1967-9",
              between(date, as.Date("1972-02-01"), as.Date("1974-07-01")) ~ "1972-74",
              between(date, as.Date("1977-01-01"), as.Date("1980-04-01")) ~ "1977-80",
              between(date, as.Date("1980-07-01"), as.Date("1981-01-01")) ~ "1980-81",
              between(date, as.Date("1983-02-01"), as.Date("1984-08-01")) ~ "1983-4", 
              between(date, as.Date("1988-03-01"), as.Date("1989-04-01")) ~ "1988-9",
              between(date, as.Date("1993-12-01"), as.Date("1995-04-01")) ~ "1993-5",
              between(date, as.Date("1999-01-01"), as.Date("2000-07-01")) ~ "1999-2000",
              between(date, as.Date("2004-05-01"), as.Date("2006-07-01")) ~ "2004-6",
              between(date, as.Date("2015-11-01"), as.Date("2019-01-01")) ~ "2015-19",
              between(date, as.Date("2022-01-01"), as.Date("2023-05-01")) ~ "2022 - ?")) %>% 
  group_by(tightening) %>% 
  mutate(small = 1 + small,
         small = (small^(12/length(small))) - 1,
         large = 1 + large,
         large = (large^(12/length(large))) - 1) %>% 
  filter(row_number() == n()) %>% 
  ungroup %>% 
  select(tightening, small, large) %>% 
  mutate(diff = small - large) %>% 
  drop_na %>% 
   kbl(digits = 3) %>% 
  kable_classic(full_width = F,
              font_size = 20,
              html_font = "serif") %>% 
              footnote(general = "\nSource: Robert Shiller data, CFA Institute",
                       footnote_as_chunk = T,
                       general_title = "")

tightening	small	large	diff
1965-6	0.039	0.008	0.032
1967-9	0.041	0.026	0.015
1972-74	-0.009	-0.030	0.022
1977-80	0.020	0.014	0.007
1980-81	0.036	-0.071	0.106
1983-4	0.062	0.068	-0.006
1988-9	0.024	0.044	-0.020
1993-5	0.025	0.021	0.004
1999-2000	-0.020	-0.010	-0.011
2004-6	-0.015	0.003	-0.018
2015-19	0.030	0.024	0.006
2022 - ?	-0.009	0.003	-0.012
Source: Robert Shiller data, CFA Institute

Financial conditions

nfci <- 
  c("NFCI") %>% 
  tq_get(get = "economic.data", from = "1920-01-01") %>% 
  select(date, nfci = price)

df %>% 
  right_join(nfci, by = "date") %>% 
  mutate(nfci_chg = nfci - lag(nfci, 1)) %>% 
  lm(large ~ nfci_chg, data = .) %>% 
  summary()

## 
## Call:
## lm(formula = large ~ nfci_chg, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.100658 -0.019118 -0.000621  0.023402  0.102335 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.013657   0.003719   3.673 0.000409 ***
## nfci_chg    -0.000693   0.004588  -0.151 0.880298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03545 on 89 degrees of freedom
##   (2703 observations deleted due to missingness)
## Multiple R-squared:  0.0002562,  Adjusted R-squared:  -0.01098 
## F-statistic: 0.02281 on 1 and 89 DF,  p-value: 0.8803

CFAI blog: Stocks and rates

Load libraries and personal ggplot theme

Load data

Correlations (stocks ~ rates and lag(rates))

Chart 1 left and right

Table 1 and Table 2(comment in filter for Table 2): by decade

Chart 2: Rolling correlations

Chart 3: Partial correlations

Table 3: Fed tightenings

Financial conditions