PREDICT 413 - Discussion Week 6

The goal for this week is to forecast the right hand side variables in a time series model using the known response variable.

The dataset consists of share prices for Safeway Inc. and Google Trends data for the word “recipe” over the past two years. Monthly share prices are used to match up with the monthly format of the Google Trends data. Because of the way Google Trends data is calculated (self-indexed within the time period and search term set), the search window is all years of Google Trends (2004-present) and only data for the term “recipe” is returned.

library(quantmod)
library(forecast)
library(lubridate)
getSymbols('KR', src = 'yahoo', 
           from = Sys.Date() - years(5), to = Sys.Date())

KR <- data.frame(
  KR,
  date = as.Date(rownames(data.frame(KR)))
  )

library(gtrendsR)
trends.dat <- gtrends('recipe', 'US')$interest_over_time

kr.rounded <- KR%>%
    left_join(
      aggregate(KR["date"], list(month = substr(KR$date, 1, 7)), max)%>%
        mutate(keep = 1)
    )%>%
  filter(keep == 1)%>%
  mutate(month = as.Date(paste0(month,'-01')))

gt.rounded <- trends.dat%>%
  mutate(month = substr(as.character(date),1,7))%>%
  group_by(month)%>%
  summarise(hits = mean(hits))%>%
  mutate(month = as.Date(paste0(month,'-01')))

kr <- inner_join(
  kr.rounded, gt.rounded,
  by = 'month'
)%>%
  select(kr.close = KR.Adjusted, month, hits)%>%
  ts(frequency = 12)

As a first exploratory step, I want to plot the decomposition of each time series to get a general idea of which components may correlate. Kroger share prices show seasonal spikes around winter and troughs in the summer and within the time period there is a generally positive trend (although it may be turning downward!). Search results for “recipe” show a downward trend over the past five years with similar but more severe seasonality around winter.

library(gridExtra)
library(ggplot2)
# devtools::install_github('yaztheme','joshyazman')
library(yaztheme)
dec.hits <- decompose(kr[,'hits'])
dec.kr <- decompose(kr[,'kr.close'])
grid.arrange(
  autoplot(dec.hits)+
    labs(title = 'Relative Search Frequency: "Recipe"',
         x = 'Years Ago', y = element_blank())+
    theme_yaz()+
    scale_x_continuous(breaks = c(1,2,3,4,5), labels = c('5','4','3','2','1'))+
    scale_color_manual(values = yaz_cols),
  autoplot(dec.kr)+
    labs(title = 'Adjusted Closing Price: "Kroger"',
         x = 'Years Ago', y = element_blank())+
    theme_yaz()+
    scale_x_continuous(breaks = c(1,2,3,4,5), labels = c('5','4','3','2','1'))+
    scale_color_manual(values = yaz_cols),
  nrow = 1
)

Using the trend in Kroger share prices to forecast average monthly search interest for “recipes” appears to require a quadratic trend model due to the curved shape of the trend line. Another option is to use the logarithm to linearize share price trends.

The first model is the quadratic combination of seasonally adjusted share prices and seasonally adjusted average search interest. This model has a root-mean-squared-error of 4.45 (for data indexed to a maximum value of 100) with residuals close to mean zero. The Breusch-Godfrey test returned a p-value of .001 indicating a low likelihood of serial autocorrelation.

quad.lm <- tslm(formula = dec.hits$trend~dec.kr$trend + dec.kr$trend^2)
quad.acc <- data.frame(accuracy(quad.lm))%>%mutate(method = 'Quadratic')

longer object length is not a multiple of shorter object length

checkresiduals(quad.lm)

    Breusch-Godfrey test for serial correlation of order up to 24

data:  Residuals from Linear regression model
LM test = 48.906, df = 24, p-value = 0.001946

The next model attempts to log-linearize seasonally adjusted share prices for Kroger and use that variable as a regressor against search interest. In this case, RMSE ticked up a notch as did the Breusch-Godfrey p-value.

fit.lm <- tslm(formula = dec.hits$trend~log(dec.kr$trend))
fit.acc <- data.frame(accuracy(fit.lm))%>%mutate(method = 'Simple')

longer object length is not a multiple of shorter object length

checkresiduals(fit.lm)

    Breusch-Godfrey test for serial correlation of order up to 24

data:  Residuals from Linear regression model
LM test = 48.372, df = 24, p-value = 0.002269

A final comparison of the two above models illustrates that the two models are very, very close in performance. That said, across all five bias and variance metrics presented below, the quadratic model ekes out a win.

mod.health <- bind_rows(quad.acc, fit.acc)%>%select(-MASE)
library(reshape2)
mod.health%>%
  melt(id.vars = 'method')%>%
  ggplot(aes(x = variable, y = value, fill = method))+
    geom_bar(stat = 'identity', position = 'dodge')+
    scale_fill_manual(name = 'Method', values = yaz_cols[3:4])+
    theme_yaz()+
    labs(title = 'Model Accuracy Measures',
         x = element_blank(), y = 'Measure Value')