R markdown for Presentation slides

Data wrangling

library(tidyquant)
library(dma)
library(tibbletime)
library(boot)
library(viridis)
library(cowplot)

# Get data via St Louis Federal Reserve's FRED database
tickers<- c("USREC",  # NBER Recession indicator
            "GS10",   # 10-year constant maturity Treasury yield
            "TB3MS",  # 3-month Tbill rate
            "PAYEMS" # nonfarm payroll employment
            )

df <- tq_get(tickers,get="economic.data",from="1945-01-01")

#variables that need transformation
my_trans <- function(in.data,transform="pctdiff3"){
  switch(transform,
         logdiff  = c(NA,diff(log(in.data))),
         pctdiff3 = 100*Delt(in.data,k=3),     # only use pctdiff3 here, but could try others
         logdiff3 = c(rep(NA,3),diff(log(in.data),3))
         )}

# variables transformed to 3  month diff:
vlist<- c("PAYEMS")

df2<- df %>% group_by(symbol) %>%
  mutate(x=ifelse(symbol %in% vlist, my_trans(price),price))

df2 %>% select(-price)  %>% spread(symbol,x) %>%
  mutate(SLOPE=GS10-TB3MS,
         CURVE=2*GS2-TB3MS-GS10,
         REC12=lead(USREC,12)) -> df3
df4 <- filter(df3, year(date)>1955 & year(date)<2016)

# recession df (for plotting)
recessions.df = read.table(textConnection(
  "Peak, Trough
  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)
}

ggplot(data=df4, aes(x=date,y=PAYEMS))+
geom_rect(data=recessions.df, inherit.aes=FALSE,
          aes(xmin=Peak, xmax=Trough, ymin=-Inf, ymax=+Inf), 
          fill='lightblue', alpha=0.25)+theme_tq()+geom_line(color="royalblue")

# fit binary model (contemporaneous)
df5 <- filter(df3, year(date)>1955 & !is.na(USREC))

xvar <- as.matrix(df5 %>% select(SLOPE,PAYEMS))
yvar <- as.matrix(df5$USREC)

mmat<- matrix(c(1,0,
                0,1,
                1,1),3,2,byrow=TRUE)

dma.fit0<- logistic.dma(unname(xvar), yvar, mmat, lambda=0.99, alpha=0.99,
                       autotune=TRUE, initialsamp=120)
# append model results 
df50 <- df5 %>% mutate(yhat0=dma.fit0$yhatdma)

# function to get coefficients out
myf<-function(j=1,                                     # model we extract
              indf=dma.fit0,                           # dma output list
              in.date=df5$date,                        # dates 
              varnames=c("INTERCEPT","SLOPE","PAYEMS") # list of variable names
){
  out=data.frame(indf$theta[j,,]) # convert array to data frame
  colnames(out)<-varnames         # rename columns
  out$date<-in.date              # dates 
  out$model<-paste0("Model:",j)   # add model label
  return(out)       
  }
# function to get cofficient variances
myfv<-function(j=1,                                     # model we extract
               indf=dma.fit0,                           # dma output list
               varnames=c("INTERCEPTv","SLOPEv","PAYEMSv") # list of variable names (v inciates variance)
               
){
  out=data.frame(indf$vartheta[j,,]) # convert array to data frame
  colnames(out)<-varnames            # rename columns
  out$model<-paste0("Model:",j)     # add model label
  return(out)
}

#use purrr::map to stack up results
data.frame(modeln=1:nrow(mmat)) %>% 
  mutate(theta  = map(modeln,myf),
         thetav = map(modeln,myfv)) %>%
  unnest(theta,thetav) %>% 
  select(model,date,c("INTERCEPT","SLOPE","PAYEMS"),
         c("INTERCEPTv","SLOPEv","PAYEMSv")) -> df.cv

# get coefficients

df.cv %>% select(model,date,c("INTERCEPT","SLOPE","PAYEMS")) %>%
  gather(var,est,-date,-model) -> df.coef

# get variances
df.cv %>% select(model,date,c("INTERCEPTv","SLOPEv","PAYEMSv")) %>%
  gather(var,v,-date,-model) %>%
  mutate(var=substr(var,1,nchar(var)-1)) -> df.var

# merge back together
df.plot<-merge(df.coef,df.var,by=c("model","date","var"))

df5 <- filter(df3, year(date)>1955 & !is.na(REC12))

xvar <- as.matrix(df5 %>% select(SLOPE,PAYEMS))
yvar <- as.matrix(df5$REC12)

mmat<- matrix(c(1,0,
                0,1,
                1,1),3,2,byrow=TRUE)

dma.fit12<- logistic.dma(unname(xvar), yvar, mmat, lambda=0.99, alpha=0.99,
                         autotune=TRUE, initialsamp=120)

# append model forecasts
df5 <- df5 %>% mutate(yhat12=dma.fit12$yhatdma)
df6<-full_join(df5,df50 %>% select(date,yhat0), by="date")
df12 <- filter(df3, year(date)>1955 & !is.na(REC12))
data.frame(modeln=1:nrow(mmat)) %>% 
  mutate(theta  = map(modeln,myf,indf=dma.fit12,in.date=df12$date),
         thetav = map(modeln,myfv,indf=dma.fit12)) %>%
  unnest(theta,thetav) %>% 
  select(model,date,c("INTERCEPT","SLOPE","PAYEMS"),
         c("INTERCEPTv","SLOPEv","PAYEMSv")) -> df.cv12

# get coefficients
df.cv12 %>% select(model,date,c("INTERCEPT","SLOPE","PAYEMS")) %>%
  gather(var,est,-date,-model) -> df.coef12

# get variances
df.cv12 %>% select(model,date,c("INTERCEPTv","SLOPEv","PAYEMSv")) %>%
  gather(var,v,-date,-model) %>%
  mutate(var=substr(var,1,nchar(var)-1)) -> df.var12

# merge back together
df.plot12<-merge(df.coef12,df.var12,by=c("model","date","var"))


# save data, called later in the .Rmd file
save(df.plot,df.plot12,dma.fit0,dma.fit12,df6,df5,df4,df3,df12,recessions.df,file="data/dma_oct2017.Rdata")