r_final_busf.R

options(digits = 3, scipen = 9999, stringasFactors = FALSE, knitr.table.format = "rst" )
remove(list = ls())
graphics.off()
#load libraries
suppressWarnings({
  suppressPackageStartupMessages({
    library(tidyverse)
    library(pdfetch)
    library(lubridate)
    library(ggplot2)
    library(forecast)
    #library(tsbox)
    library(tsibble)
    #library(timetk)
    library(TSstudio)
    
    library(rio)
    library(pdfetch)
    
    library(tidyr)
    library(stringr)
    library(dygraphs)
    library(quantmod)
    
    library(broom)
  })
})

rwct=pdfetch::pdfetch_FRED("RTWVDCT684NMFRBDAL")
emp=pdfetch::pdfetch_FRED("CTNA")
earn=pdfetch::pdfetch_FRED("SMU09000000500000003")
hours=pdfetch::pdfetch_FRED("SMU09000000500000002")

ts_info(rwct)

##  The rwct series is a xts object with 1 variable and 426 observations
##  Frequency: monthly 
##  Start time: 1988-01-01 
##  End time: 2023-06-01

ts_info(emp)

##  The emp series is a xts object with 1 variable and 411 observations
##  Frequency: monthly 
##  Start time: 1990-01-01 
##  End time: 2024-03-01

ts_info(earn)

##  The earn series is a xts object with 1 variable and 207 observations
##  Frequency: monthly 
##  Start time: 2007-01-01 
##  End time: 2024-03-01

ts_info(hours)

##  The hours series is a xts object with 1 variable and 207 observations
##  Frequency: monthly 
##  Start time: 2007-01-01 
##  End time: 2024-03-01

rwct.ts=ts(rwct,start=c(1988,1),frequency = 12)
rwct.ts=window(rwct.ts,start=c(2007,1),end=c(2023,6),frequency=12)
emp.ts=ts(emp,start=c(1990,1),frequency=12)
emp.ts=window(emp.ts,start=c(2007,1),end=c(2023,6),frequency=12)

earn.ts=ts(earn,start=c(2007,1),frequency=12)
earn.ts=window(earn.ts,start=c(2007,1),end=c(2023,6),frequency=12)
hours.ts=ts(hours,start=c(2007,1),frequency=12)
hours.ts=window(hours.ts,start=c(2007,1),end=c(2023,6),frequency=12)

cbind(rwct.ts,emp.ts,earn.ts,hours.ts)|>tidy()|>
  group_by(series) |>
  dplyr::mutate(normindex=value*100/first(value),
                standvalues=value/max(value))|>
  ggplot(aes(x=index,y=normindex,col=series))+geom_line()
#ggplot(aes(x=index,y=standvalues,col=series))+geom_line()

#standarize the variable 
mydat=cbind.data.frame(rwct.ts/max(rwct.ts),
                       emp.ts/max(emp.ts),
                       earn.ts/max(earn.ts),
                       hours.ts/max(hours.ts)
)
mydat=mydat|>dplyr::mutate(CTIndex=rowSums(mydat))
head(mydat,3)

##   RTWVDCT684NMFRBDAL  CTNA SMU09000000500000003 SMU09000000500000002 CTIndex
## 1              0.780 0.987                0.718                0.980    3.47
## 2              0.778 0.986                0.729                0.968    3.46
## 3              0.777 0.986                0.722                0.977    3.46

mydat.ts=ts(mydat$CTIndex,start=c(2007,1),frequency = 12)

mydat.ts|>tidy()|>
  ggplot(aes(x=index,y=value))+geom_line()

#forecast #######

mymodel=HoltWinters(mydat.ts)
myfc=forecast(mymodel,h=24)
autoplot(myfc)
autoplot(mydat.ts)+autolayer(myfc$mean) + 
  xlim(2022,2026)+ylim(3.7,4.0)

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.