Libraries and theme and NBER recession dates
# Libraries
library(tidyquant)
library(tidyverse)
library(broom)
library(knitr)
library(psych)
library(data.table)
library(janitor)
library(magrittr)
library(sandwich)
library(car)
library(slider)
library(lmtest)
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"))
recessions_df =
read.table(
textConnection(
"begin, end
1857-06-01, 1858-12-01
1860-10-01, 1861-06-01
1865-04-01, 1867-12-01
1869-06-01, 1870-12-01
1873-10-01, 1879-03-01
1882-03-01, 1885-05-01
1887-03-01, 1888-04-01
1890-07-01, 1891-05-01
1893-01-01, 1894-06-01
1895-12-01, 1897-06-01
1899-06-01, 1900-12-01
1902-09-01, 1904-08-01
1907-05-01, 1908-06-01
1910-01-01, 1912-01-01
1913-01-01, 1914-12-01
1918-08-01, 1919-03-01
1920-01-01, 1921-07-01
1923-05-01, 1924-07-01
1926-10-01, 1927-11-01
1929-08-01, 1933-03-01
1937-05-01, 1938-06-01
1945-02-01, 1945-10-01
1948-11-01, 1949-10-01
1953-07-01, 1954-05-01
1957-08-01, 1958-04-01
1960-04-01, 1961-02-01
1969-12-01, 1970-11-01
1973-11-01, 1975-03-01
1980-01-01, 1980-07-01
1981-07-01, 1982-11-01
1990-07-01, 1991-03-01
2001-03-01, 2001-11-01
2007-12-01, 2009-06-01
2020-04-01, 2020-06-01"
),
sep = ',',
colClasses = c('Date', 'Date'),
header = TRUE
)
Ex. 1: real spot price of gold
cpi <-
c("CPIAUCSL") %>%
tq_get(get = "economic.data", from = "1978-12-01") %>%
select(date, cpi = price)
gold <-
read_csv("gold.csv") %>%
mutate(date = date(mdy_hm(date)),
date = floor_date(date, "month")) %>%
group_by(date) %>%
filter(row_number() == n()) %>%
ungroup
gold %>%
left_join(cpi, by = "date") %>%
mutate(real_gold = 100*(gold/cpi)) %>%
drop_na %>%
ggplot(aes(x = date, y = real_gold)) +
geom_line() +
geom_text(aes(label = paste0("$", round(real_gold, 1))), data = data.table::last, vjust = -1) +
geom_point(data = data.table::last, color = "gold", size = 4) +
theme +
labs(title = "Real spot price of gold",
caption = "Note: Chart shows spot price of gold deflated using the Consumer Price Index level.\nSource: World Gold Council, FRED",
subtitle = "1979 - March 2024",
x = NULL,
y = "Real gold spot price")
Dataframe for other exhibits
median <-
read_csv("median_may.csv") %>%
mutate(date = mdy(date),
median = median_pce/100) %>%
select(date, median)
gold <-
read_csv("gold.csv") %>%
mutate(date = date(mdy_hm(date)),
date = floor_date(date, "month")) %>%
group_by(date) %>%
filter(row_number() == n()) %>%
ungroup %>%
mutate(gold = (gold/lag(gold, 1)) - 1)
core_pce <-
c("PCEPILFE") %>%
tq_get(get = "economic.data", from = "1978-12-01") %>%
select(date, core_pce = price) %>%
mutate(core_pce = (core_pce/lag(core_pce, 1)) - 1)
pce <-
c("PCEPI") %>%
tq_get(get = "economic.data", from = "1978-12-01") %>%
select(date, pce = price) %>%
mutate(headline_pce = (pce/lag(pce, 1)) - 1) %>%
select(date, headline_pce)
cpi <-
c("CPIAUCSL") %>%
tq_get(get = "economic.data", from = "1978-12-01") %>%
select(date, headline_cpi = price) %>%
mutate(headline_cpi = (headline_cpi/lag(headline_cpi, 1)) - 1) %>%
select(date, headline_cpi)
core_cpi <-
c("CPIAUCSL") %>%
tq_get(get = "economic.data", from = "1978-12-01") %>%
select(date, core_cpi = price) %>%
mutate(core_cpi = (core_cpi/lag(core_cpi, 1)) - 1) %>%
select(date, core_cpi)
df <-
gold %>%
left_join(median, by = "date") %>%
left_join(core_pce, by = "date") %>%
left_join(pce, by = "date") %>%
left_join(core_cpi, by = "date") %>%
left_join(cpi, by = "date")
drop_na
## function (data, ...)
## {
## check_dots_unnamed()
## UseMethod("drop_na")
## }
## <bytecode: 0x00000272c8ad03a0>
## <environment: namespace:tidyr>
Corr and regression (gold ~ headline_pce)
df %>%
#cor.test(~headline_pce + gold, data = .)
lm(gold ~ headline_pce, data = .) %>% summary
##
## Call:
## lm(formula = gold ~ headline_pce, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.246275 -0.029965 -0.001953 0.026138 0.255919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00149 0.00315 0.473 0.6364
## headline_pce 1.71595 0.92126 1.863 0.0631 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05218 on 542 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.00636, Adjusted R-squared: 0.004527
## F-statistic: 3.469 on 1 and 542 DF, p-value: 0.06306
#coeftest(vcov = vcovHAC)
Ex. 2: Scatterplot of headline_pce vs. gold
df %>%
ggplot(aes(x = headline_pce, y = gold)) +
geom_point(shape = 22, fill = "gold4") +
theme +
geom_smooth(method = "lm", se = T) +
labs(x = "Headline PCE deflator m/m % change", y = "Gold spot price m/m % change",
title = "Headline PCE vs. gold",
subtitle = "Monthly changes, 1979 to March 2024",
caption = "Note: Chart shows pairs of monthly changes in headline PCE inflation and gold spot price with linear regression best fit line\n and shaded standard error values.\nSource: FRED, World Gold Council")
Robustness checks
df %>%
select(date, gold, headline_cpi) %>%
cor.test(~ gold + headline_cpi, data = .)
##
## Pearson's product-moment correlation
##
## data: gold and headline_cpi
## t = 2.3598, df = 542, p-value = 0.01864
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01692268 0.18335785
## sample estimates:
## cor
## 0.1008458
#lm(gold ~ headline_cpi, data = .) %>%
summary
## standardGeneric for "summary" defined from package "base"
##
## function (object, ...)
## standardGeneric("summary")
## <environment: 0x00000272ca44f160>
## Methods may be defined for arguments: object
## Use showMethods(summary) for currently available ones.
Ex. 3: Rolling regression - headline inflation
df <-
df %>%
mutate(index = row_number())
# assign window length in index units (months in this example)
window <- 36
# initialize vector of length of periods over which regressions will be run (= total length minus window length)
results_vector <-
tibble(values = rep(0, length(df$index) - window),
errors = rep(0, length(df$index) - window)) # getting values for total less observations lost due to window length
for(i in 1:length(results_vector$values)){
mod <-
df %>%
filter(between(index, first(index) + i - 1, first(index) + window + i - 1)) %>%
lm(gold ~ headline_pce, data = .) %>%
coeftest
results_vector[i,1] <- mod[2, 1] #extract coefficient
results_vector[i, 2] <- mod[2, 2] #extract SE
}
results_vector %>%
mutate(date = df$date[df$index > window]) %>% # Add date column by subsetting original date vector
ggplot() +
geom_ribbon(aes(x = date, ymin = values - 2*errors, ymax = values + 2*errors), fill = "gold") +
geom_line(aes(x = date, y = values)) +
geom_rect(aes(xmin = begin, xmax = end, ymin = -Inf, ymax = +Inf),
fill = "azure3", data = recessions_df %>% filter(year(begin) > 1978), alpha = .5) +
geom_hline(yintercept = 0, linetype = 2) +
theme +
labs(title = "Gold 36-month rolling headline-inflation beta",
subtitle = "1979 - March 2024",
caption = "Note: Chart shows the coefficient estimated from a regression of monthly change in gold spot price\non monthly change on median PCE inflation with 2-standard-error band.\n NBER recession dates shaded in gray.\nSource: World Gold Council, FRED, NBER",
x = NULL,
y = "Gold spot price inflation beta")
lm(results_vector$values ~ lag(results_vector$values, 1)) %>% summary
##
## Call:
## lm(formula = results_vector$values ~ lag(results_vector$values,
## 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5467 -0.5846 -0.0406 0.4921 7.7910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07652 0.07044 1.086 0.278
## lag(results_vector$values, 1) 0.94701 0.01431 66.175 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.497 on 507 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8962, Adjusted R-squared: 0.896
## F-statistic: 4379 on 1 and 507 DF, p-value: < 2.2e-16
Other test: gold and economic slack
gold <-
read_csv("gold.csv") %>%
mutate(date = date(mdy_hm(date)),
date = floor_date(date, "quarter")) %>%
group_by(date) %>%
filter(row_number() == n()) %>%
ungroup
rgdp <-
c("GDPC1") %>%
tq_get(get = "economic.data", from = "1978-01-01") %>%
select(date, gdp = price)
potgdp <-
c("GDPPOT") %>%
tq_get(get = "economic.data", from = "1978-01-01") %>%
select(date, gdppot = price)
gold %>%
left_join(rgdp, by = "date") %>%
left_join(potgdp, by = "date") %>%
mutate(gap = gdp - gdppot,
gap = (gap / lag(gap, 1)) - 1,
gold_chg = (gold / lag(gold, 1)) - 1) %>%
lm(gold_chg ~ gap, data = .) %>%
summary
##
## Call:
## lm(formula = gold_chg ~ gap, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27082 -0.05213 -0.00889 0.04904 0.41525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0165524 0.0063788 2.595 0.0102 *
## gap 0.0010275 0.0009705 1.059 0.2911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08568 on 179 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.006224, Adjusted R-squared: 0.0006718
## F-statistic: 1.121 on 1 and 179 DF, p-value: 0.2911
