Using a time series model to inspect the seasonally adj. US CPI inflation measure, % change

inflation.R

library(tidyverse)
library(forecast)
library(urca)
library(dynlm)
library(lmtest)
library(sandwich)

# https://data.bls.gov/timeseries/CUSR0000SA0&output_view=pct_1mth
read_csv("~/Downloads/inflation.csv", skip = 11) %>%
  janitor::clean_names() %>%
  mutate(period = gsub("M", "", period),
         date = paste(year, period, "01", sep = "-") %>% as.POSIXct(),
         cpi = value) -> d

acf(d$cpi)
pacf(d$cpi)

auto.arima(
  d$cpi %>% tail(120), # last 10 yrs to focus on recent dynamics
  seasonal = FALSE, # already seasonally adj
  approximation = FALSE
) -> fit

forecast(fit, 1) %>% plot(main = "1 month forecast")
abline(h = 0, col = "red")
plot(resid(fit), type = "b",
     pch = 20,
     main = "Discrepancy between actuals and the model")
abline(h = 0, col = "red")

# https://fred.stlouisfed.org/series/MTSO133FMS
read_csv("~/Downloads/MTSO133FMS.csv") %>%
  janitor::clean_names() %>%
  rename(outlays = mtso133fms) %>%
  mutate(log_outlays = log(outlays)) %>%
  mutate(diff_log_outlays = c(NA, diff(log_outlays))) %>%
  filter(!is.na(diff_log_outlays)) -> outlays

fracdiff::fracdiff(d$cpi) %>% summary()
tseries::kpss.test(d$cpi)
auto.arima(d$cpi)
outlays %>%
  inner_join(d) %>%
  left_join(
    read_csv("~/Downloads/data_table_for_daily_case_trends__the_united_states.csv", skip = 2) %>%
      janitor::clean_names() %>%
      mutate(date = as.POSIXct(date, format = "%b %d %Y")) %>%
      group_by(date = lubridate::floor_date(date, "month")) %>%
      # we have averaged about 3.2M cases per month
      summarise(covid_cases = sum(new_cases)/1e6/3.2) %>%
      select(date, covid_cases) %>%
      replace_na(list(covid_cases = 0)) %>%
      arrange(date)
  ) %>%
  replace_na(list(covid_cases = 0)) %>%
  mutate(outlays = outlays / 1e6,
         log_outlays = log(outlays),
         trend = 1:n()) -> d2

ca.jo(d2 %>% select(cpi, log_outlays, covid_cases), ecdet = "none") %>% summary()
MASS::rlm(cpi ~ lag(log_outlays) + lag(covid_cases), data = d2) -> fit
fracdiff::fracdiff(resid(fit))
summary(fit)
par(mfrow = c(1, 1)); plot(resid(fit), type = "b"); abline(h = 0, col = "red")

lm(cpi ~ log_outlays_ts + covid_cases) -> fit0
summary(fit0)
coeftest(fit0, NeweyWest(fit0))


library(dynlm)
cpi <- ts(d2$cpi)
log_outlays_ts <- ts(d2$log_outlays)
covid_cases <- ts(d2$covid_cases)
dynlm(d(cpi) ~ L(cpi) + L(d(cpi), 1:3) +
        L(log_outlays_ts) + L(d(log_outlays_ts), 1:12) +
        L(covid_cases) + L(d(covid_cases), 1:6) - 1) -> fit
# acf(resid(fit))
bgtest(fit)
bptest(fit)
fracdiff::fracdiff(resid(fit))
summary(fit)
auto.arima(resid(fit), seasonal = F)
par(mfrow = c(1, 1)); plot(ts(resid(fit), frequency = 12, start = c(1980, 11)), type = "l"); abline(h = 0, col = "red")
coeftest(fit, NeweyWest(fit))


read_csv("~/Downloads/DGS1.csv") %>%
  janitor::clean_names() %>%
  inner_join(d2) %>%
  mutate(diff_log_dgs1 = c(NA, diff(log(dgs1))),
         log_dgs1 = log(dgs1)) %>%
  filter(!is.na(diff_log_dgs1)) -> d3

plot(d3$dgs1, type = "l")
plot(diff(log(d3$dgs1)), type = "l")

cpi <- ts(d3$cpi)
log_outlays_ts <- ts(d3$log_outlays)
covid_cases <- ts(d3$covid_cases)
log_one_yr_treasuries_interest_rate <- ts(d3$dgs1)
auto.arima(log_one_yr_treasuries_interest_rate)

dynlm(d(cpi) ~ L(cpi) + L(d(cpi)) +
        L(log_outlays_ts) + L(d(log_outlays_ts), 1:12) +
        L(covid_cases) + L(d(covid_cases), 1:6) - 1) -> fit0
dynlm(d(cpi) ~ L(cpi) + L(d(cpi)) +
        L(log_outlays_ts) + L(d(log_outlays_ts), 1:12) +
        L(covid_cases) + L(d(covid_cases), 1:12) - 1) -> fit00
BIC(fit00, fit0)
bgtest(fit0)
bptest(fit0)
Box.test(resid(fit0), type = "Ljung-Box")
tseries::kpss.test(resid(fit0))
fracdiff::fracdiff(resid(fit0))
summary(fit0)
auto.arima(resid(fit0), seasonal = F)
acf(resid(fit0))
par(mfrow = c(1, 1)); plot(ts(resid(fit0), frequency = 12, start = c(1980, 11)), type = "l"); abline(h = 0, col = "red")
coeftest(fit0, NeweyWest(fit0))