First load required packages.

library(tidyverse)
library(lubridate)
library(wavelets)
library(forecast)
library(tidyquant)
library(data.table)
library(knitr)
library(kableExtra)
# not needed but could be useful
#library(wmtsa)

Our data will come in three flavors. First we’ll create a small dataframe of NBER business cycle dates by hand, then we’ll grab some monthly time series and finally we’ll grab some quarterly time series. The wavelets are a bit more intuitive at quarterly frequency, so we’ll create quarterly averages of the monthly data. Note that this will create some slight differences between our results and Matthes, Lubik and Verona (2019) because they aggregated daily data into quarterly data for the interest rate series. Subsequent data revisions and different time periods should explain most of the other (very small) differences in what we’re about to produce.

Data pulled on April 11, 2019

# recession data ----

recessions.df = read.table(textConnection(
  "Peak, Trough
  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"), sep=',',
  colClasses=c('Date', 'Date'), header=TRUE)%>%
  # these dates at numbers will be helpful with ts() plots laters
  mutate(Peakt = year(Peak)+(month(Peak)-1)/12,
         Trought = year(Trough)+(month(Trough)-1)/12)




# fred data ----

# monthly data----

# FRED tickers
tickers=c("UNRATE","FEDFUNDS","GS10",
          "HSN1F","MSACSR",
          "HOUST","HOUST1F",
          "HOUST5F", "TTLHHM156N") 

# Human readable descriptions
lookup_1 <- c("Unemployment Rate (%)",
              "Federal Funds Effective Rate (%)",
              "10-Year Treasury Constant Maturity Rate (%)",
              "New One Family Houses Sold (1000s SAAR)",
              "Monthly Supply of New Houses for Sale (Months)",
              "Housing Starts (total, 1000s SAAR)",
              "Housing Starts (1-unit structures, 1000s SAAR)",
              "Housing Starts: (5-Unit Structures or More, 1000s SAAR)", 
              "Household Estimates (Thousands)"
              )


dfm <- tidyquant::tq_get(tickers,get="economic.data",from="1954-01-01")

# convert monthly data to quarter data
dfq <- 
  dfm %>% mutate(year=year(date), quarter=(month(date)-1) %/% 3 + 1) %>%
  group_by(symbol,year,quarter) %>%
  summarize(vq=mean(price),
            date=first(date)) %>%
  ungroup()

# spread it out (wide is more convenient than tidy below)

dfq2 <- dfq %>% select(symbol,vq,date) %>% group_by(date) %>% spread(symbol,vq)

# quarterly data ----

# FRED tickers

tickers2 <- c("RRVRUSQ156N",
              "RHVRUSQ156N",
              "USSTHPI",
              "PRFI",  
              "A011RL1Q225SBEA",
              "A011RY2Q224SBEA"
       
)

# Human readable descriptions
lookup_2 <- c( 
  "Rental Vacancy Rate for the United States (%)",
  "Homeowner Vacancy Rate for the United States (%) ",
  "FHFA All-Transactions House Price Index for the United States (Index 1980 Q1 = 100)",
  "Nominal private residential fixed investment ($ Billions, SAAR)",
  "Percent change in real private residential fixed investment (%)",
  "Contribution to GDP growth real private residential fixed investment (%)")

# lookup table (tickers and series_name)
df_lookup=data.frame(series=c(tickers,tickers2), series_name=c(lookup_1,lookup_2), stringsAsFactors = FALSE)

# quarterly
# spread it out (wide is more convenient than tidy below)
dfq3<- tidyquant::tq_get(tickers2,get="economic.data",from="1947-01-01") %>% group_by(date) %>% spread(symbol,price)


dfq4 <- merge(dfq2, dfq3, all="TRUE", by="date") %>% mutate(tt=year(date)+ (month(date)-1)/12) %>% ungroup()

# create quarerly growth in households, house prices

dfq4 <- mutate(dfq4, 
               HPA= log(USSTHPI)-lag(log(USSTHPI)),
               HHG = log(TTLHHM156N)- lag(log(TTLHHM156N))
               )

# adjust lookup for new transformed variables
df_lookup = rbind(df_lookup, data.frame(series=c("HPA","HHG"),
                                        series_name=c("Quarterly House Price Growth (%)","Quarterly Household Growth (%)"),
                                        stringsAsFactors = FALSE))

f_ts <- function(invar="GS10"){
  invar <- enquo(invar)
  df <- select(dfq4, tt,!!invar)
  df <- df[complete.cases(df),]
  yts <- ts(df[,2], start=min(df$tt), frequency=4)
}

g <- 
autoplot(f_ts("HOUST"))+
  labs(x="date (quarterly)",y="", title="U.S. Housing Starts (1000s, SAAR)",
       caption="@lenkiefer Source: U.S. Census Bureau and U.S. Department of Housing and Urban Development, Housing Starts: Total: New Privately Owned Housing Units Started [HOUST],\nretrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/HOUST, April 11, 2019.") +
  theme(plot.caption=element_text(hjust=0))

# add recession shading

gbox <- 
  geom_rect(data=filter(recessions.df,Peakt>1954), inherit.aes=F, 
            aes(xmin=Peakt, xmax=Trought, ymin=-Inf, ymax=+Inf), fill='darkgray', alpha=0.5)

# we want to move the box in front of the autoplot, so we ce rearrange like so
# see: https://stackoverflow.com/questions/20249653/insert-layer-underneath-existing-layers-in-ggplot2-object
g$layers <- c(gbox, g$layers)

There are many options for the wavelets, but following Matthes, Lubik and Verona (2019) we will use the Haar wavelet filter with 6 levels. That will break the variation up between 2-4 quarters, 4-8 quarters, 8-16 quarters, 16-32 quarters, 32-64 quarters, 64-128 quarters and the long-term trend. We will call the first 2 groups “Short Term” the next two “Business Cycle”,and the next 2 “Medium Term”.

# function to analyze wavelets
myf <- function(invar){
  invar <- enquo(invar)
  df <- select(dfq4, tt,!!invar)
  df <- df[complete.cases(df),]
  yts <- ts(df[,2], start=min(df$tt), frequency=4)
  w2 <- modwt(yts, filter="haar", n.levels=6, boundary="periodic", fast=TRUE)
  Wdf <- data.frame(w2@W)
  var_decomp <-  c(round(100*colSums(Wdf**2)/(nrow(Wdf)-1)/ var(yts),1), 
                   long_term=round(100*(1-sum(colSums(Wdf**2)/(nrow(Wdf)-1)/ var(yts))),1))

  Sj <- yts-rowSums(Wdf)
  tWdf <- ts(Wdf, start=min(df$tt), frequency=4)
  D12 <- ts(rowSums(Wdf[,1:2]), start=min(df$tt), frequency=4)
  D34 <- ts(rowSums(Wdf[,3:4]), start=min(df$tt), frequency=4)
  D56 <- ts(rowSums(Wdf[,5:6]), start=min(df$tt), frequency=4)
  return(list(vd=var_decomp,out_ts=ts.union(D12,D34,D56,Sj),mod=w2))
}

# function to plot retults

my_plotf <- function(invar="A011RL1Q225SBEA"){
  plot_subtitle <- df_lookup[df_lookup$series==invar,"series_name"]
  invar <- enquo(invar)
  out <- myf(!!invar)
  ss<- min(time(f_ts(!!invar)))+ 2
  plot_ts2 <- window(f_ts(!!invar), start=ss)
  plot_ts <- window(out$out_ts, start=ss)
  plot_df <- data.frame(tt=time(plot_ts), plot_ts)  %>% 
    gather(series,val,-tt) %>%
    # labels for series
    mutate( series_id = case_when(series=="D12" ~ "a) Short Term: 2-8 Quarters",
                                  series=="D34" ~ "b) Business cycle: 16-32 Quarters",
                                  series=="D56" ~ "c) Medium Term: 64-128 Quarters",
                                  T ~ "d) Long Term: >128 Quarters"
    ))
  
  plot_df2 <- data.frame(tt=time(plot_ts2), plot_ts2)
  
  ggplot(data=plot_df, aes(x=tt, y=val))+
    facet_wrap(~series_id, ncol=1, scales="free_y")+
    geom_hline(aes(yintercept=ifelse(series=="Sj",NA,0)), linetype=2)+
    geom_rect(data=filter(recessions.df,Peakt>ss-1), inherit.aes=F, 
              aes(xmin=Peakt, xmax=Trought, ymin=-Inf, ymax=+Inf), fill='darkgray', alpha=0.5) +
    geom_line(size=1.05)+
    scale_y_continuous(sec.axis=dup_axis())+
    labs(subtitle=plot_subtitle,
         title="Wavelet Decomposition",
         x="date",
         y="Variation due to fluctuations at frequency:",
         caption="@lenkiefer Based on decomposition described in Matthes, C., Lubik, T. A., & Verona, F. (2019). \n'Assessing US Aggregate Fluctuations Across Time and Frequencies' Richmond Fed Working Paper (No. 19-6).\nData retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/\nGray bars NBER Recessions")+
    theme(plot.caption=element_text(hjust=0))
}

myf2 <- function(invar1="HOUST"){
  invar2 <- enquo(invar1)
  v1 <- data.frame(t(myf(!!invar2)$vd))
  return(v1)
}

# custom labels for table

df_table <- df_lookup %>% mutate(v=map(series,myf2)) %>% unnest(v) %>% select(-series) 
colnames(df_table) <- c("Series","D1:  2-4Q", "D2:  4-8Q","D3:  8-16Q","D4:  16-31Q","D5  32-64Q","D6:  64-128Q", "S6: >128Q")
# reorganize table
df_table$id <- c(1,1,1, 2,2,3,3,3,4,5,5,6,2,2,2,6,4)
df_table <- arrange(df_table,id)
df_table$id <- NULL

kable(df_table, 
      font_size = 8, 
      row.names=FALSE, 
      caption=" Variance Decomposition for U.S. Macroeconomic and Housing Market Indicators (quarterly data)") %>%
  kable_styling("striped") %>%
  add_header_above(c(" " = 1, "Short Term" = 2, "Business Cycle" = 2, "Medium Term" = 2, "Long Term" =1 ))%>%
  pack_rows(index = c( "Macro Variables " = 3, 
                       "New Home Sales" = 2, 
                       "Housing Investment"= 3, 
                       "Housing Construction"=3, 
                       "Households"=2,
                       "Vacancy Rates"=2, 
                       "House Prices" = 2)) %>%
  footnote(general="@lenkiefer Based on decomposition described in Matthes, C., Lubik, T. A., & Verona, F. (2019). \n'Assessing US Aggregate Fluctuations Across Time and Frequencies' Richmond Fed Working Paper (No. 19-6).\nData retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/\nSeries length varies by series, maximum available quarters selected",
           number="Data Retrieved on April 11, 2019")