#The 'seasona.duration' argument tells the model how many time points each season should last for. The 'nseasons' argument tells it how many seasons are in a cycle. The total number of time points in a cycle is season.duration * nseasons.
#AddSeasonal(..., nseasons = 52, season.duration = 7). for daily data
# GeneralizedLocalLinearTrend (sorry about the nondescriptive name) assumes the "slope" component of trend is an AR1 process instead of a random walk. It is my default option if I want to forecast far into the future. Most of your time series variation seems to come from seasonality, so you might try AddLocalLevel or even AddAr instead of AddLocalLinearTrend



#Bsts offers support for multiple seasonalities. For example, if you have several weeks of hourly data then you will have an hour-of-day effect as well as a day-of-week effect. You can model these using a single seasonal effect with 168 seasons (which would allow for different hourly effects on weekends and weekdays), or you can assume additive seasonal patterns using the season.duration argument to AddSeasonal,

#ss <- AddSeasonal(ss, y, nseasons = 24)
#ss <- AddSeasonal(ss, y, nseasons = 7, season.duration = 24)
library(rio)
library(stringr)
library(tseries)
library(lubridate)
library(tidyverse)
library(timetk)
library(tidyquant)
library(scales)
library(forecast)   #  forecasting pkg
library(sweep)# Broom tidiers for forecast pk
library(timekit)
library(plotly)
library(broom)
library(tibble)
library(ggfortify)  #similar to broom,also helps autoplot of lm objects
library(gam) #Fit general additive model
library(highcharter)
library(dygraphs)
library(DT)
require(RColorBrewer)
library(viridisLite)
library(bsts)
library(readxl)
select=dplyr::select

setwd("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/time series")
#load excel file with rio
UptimeData<- rio::import("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/Uptimebycategories08-29-2017b.xlsx")

DT::datatable(UptimeData) 
# replace white spaces in column names with _

names(UptimeData)<-stringr::str_replace_all(names(UptimeData),"\\s", "_")

UptimeData%>%head()%>%knitr::kable()
Month-Yr Dur/EBRG priority_category_cd build_phase Build_Category Day_of_Week vehicle_family_ref APG_Uptime
Jan-2017 Durability P70 Comp Veh Comp Veh Sun JT 0
Jan-2017 Durability P70 Comp Veh Comp Veh Mon JT 0
Jan-2017 Durability P70 Comp Veh Comp Veh Tue JT 0
Jan-2017 Durability P70 Comp Veh Comp Veh Wed JT 0
Jan-2017 Durability P70 Comp Veh Comp Veh Thu JT 0
Jan-2017 Durability P70 Comp Veh Comp Veh Fri JT 0 
UptimeData=UptimeData%>%mutate(Month=stringr:: str_sub(`Month-Yr`, 1, 3))%>%mutate_if(is.character,factor)


UptimeData%>%head()%>%knitr::kable()
Month-Yr Dur/EBRG priority_category_cd build_phase Build_Category Day_of_Week vehicle_family_ref APG_Uptime Month
Jan-2017 Durability P70 Comp Veh Comp Veh Sun JT 0 Jan
Jan-2017 Durability P70 Comp Veh Comp Veh Mon JT 0 Jan
Jan-2017 Durability P70 Comp Veh Comp Veh Tue JT 0 Jan
Jan-2017 Durability P70 Comp Veh Comp Veh Wed JT 0 Jan
Jan-2017 Durability P70 Comp Veh Comp Veh Thu JT 0 Jan
Jan-2017 Durability P70 Comp Veh Comp Veh Fri JT 0 Jan 
Durability=UptimeData%>%filter(`Dur/EBRG`=="Durability")%>%rename(Dur=`Dur/EBRG`)

DT::datatable(Durability)
summary(Durability)
##      Month-Yr           Dur      priority_category_cd   build_phase 
##  May-2017:138   Durability:838   P100:269             VP      :372  
##  Jun-2017:121   EBRG-Trans:  0   P50 :  0             PS      :113  
##  Feb-2017:112                    P70 :487             Proto   :110  
##  Jul-2017:101                    P85 : 82             Comp Veh: 66  
##  Apr-2017:100                                         J1      : 53  
##  Jan-2017: 94                                         VP(A)   : 52  
##  (Other) :172                                         (Other) : 72  
##     Build_Category Day_of_Week vehicle_family_ref   APG_Uptime    
##  Comp Veh  : 66    Fri:125     KL     :103        Min.   :0.0000  
##  J1        : 53    Mon:121     JL     : 99        1st Qu.:0.4963  
##  Mule_Proto:161    Sat:115     LA-SRT : 83        Median :0.7599  
##  PS        :113    Sun:110     JT     : 66        Mean   :0.6764  
##  VP        :438    Thu:125     MP     : 60        3rd Qu.:0.9727  
##  NA's      :  7    Tue:121     RU     : 49        Max.   :1.0000  
##                    Wed:121     (Other):378                        
##      Month    
##  May    :138  
##  Jun    :121  
##  Feb    :112  
##  Jul    :101  
##  Apr    :100  
##  Jan    : 94  
##  (Other):172
smy=summary(Durability)%>%as.data.frame.array()

#smy=smy%>%replace_na(list(Dur = "",priority_category_cd="",Build_Category="",APG_Uptime="")) 
#replace NA variables in table
#smy=smy%>% replace_na(list(Dur = "",priority_category_cd= ""))

smy
Durability=na.omit(Durability)

Durability_ts=ts(Durability$APG_Uptime)




#load excel file with rio
CPGUptimeData<- rio::import("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/CPGUptimecategories.xlsx",sheet="CPG")

DT::datatable(CPGUptimeData) 
# replace white spaces in column names with _

names(CPGUptimeData)<-stringr::str_replace_all(names(CPGUptimeData),"\\s", "_")

CPGUptimeData%>%head()%>%knitr::kable()
Month-Yr Dur/EBRG priority_category_cd build_phase Day_of_Week vehicle_family_ref CPG_Uptime
Jan-2017 Durability P100 J1 Sun RU-PHEV 0.8333333
Jan-2017 Durability P100 J1 Mon RU-PHEV 1.0000000
Jan-2017 Durability P100 J1 Tue RU-PHEV 0.6511628
Jan-2017 Durability P100 J1 Wed RU-PHEV 0.3333333
Jan-2017 Durability P100 J1 Thu RU-PHEV 1.0000000
Jan-2017 Durability P100 J1 Fri RU-PHEV 1.0000000 
CPGUptimeData=CPGUptimeData%>%mutate(Month=stringr:: str_sub(`Month-Yr`, 1, 3))%>%mutate_if(is.character,factor)


CPGUptimeData%>%head()%>%knitr::kable()
Month-Yr Dur/EBRG priority_category_cd build_phase Day_of_Week vehicle_family_ref CPG_Uptime Month
Jan-2017 Durability P100 J1 Sun RU-PHEV 0.8333333 Jan
Jan-2017 Durability P100 J1 Mon RU-PHEV 1.0000000 Jan
Jan-2017 Durability P100 J1 Tue RU-PHEV 0.6511628 Jan
Jan-2017 Durability P100 J1 Wed RU-PHEV 0.3333333 Jan
Jan-2017 Durability P100 J1 Thu RU-PHEV 1.0000000 Jan
Jan-2017 Durability P100 J1 Fri RU-PHEV 1.0000000 Jan 
CPGDurability=CPGUptimeData%>%filter(`Dur/EBRG`=="Durability")%>%rename(Dur=`Dur/EBRG`)

DT::datatable(CPGDurability)
summary(CPGDurability)
##      Month-Yr            Dur       priority_category_cd   build_phase 
##  Mar-2017:181   Durability :1122   P100:327             VP      :522  
##  Jul-2017:153   EBRG-Trans :   0   P70 :663             Proto   :264  
##  Jun-2017:151   Reliability:   0   P85 :132             PS      :135  
##  May-2017:144                                           Mule    :110  
##  Feb-2017:141                                           J1      : 30  
##  Apr-2017:124                                           Comp Veh: 21  
##  (Other) :228                                           (Other) : 40  
##  Day_of_Week vehicle_family_ref   CPG_Uptime         Month    
##  Fri:160     JL     :212        Min.   :0.0000   Mar    :181  
##  Mon:163     KL     :130        1st Qu.:0.5288   Jul    :153  
##  Sat:152     M1     :103        Median :0.7770   Jun    :151  
##  Sun:152     DT     : 92        Mean   :0.7023   May    :144  
##  Thu:163     MP     : 83        3rd Qu.:0.9725   Feb    :141  
##  Tue:166     RU-PHEV: 80        Max.   :1.0000   Apr    :124  
##  Wed:166     (Other):422                         (Other):228
smy=summary(CPGDurability)%>%as.data.frame.array()

#smy=smy%>%replace_na(list(Dur = "",priority_category_cd="",Build_Category="",APG_Uptime="")) 
#replace NA variables in table
#smy=smy%>% replace_na(list(Dur = "",priority_category_cd= ""))

smy%>%knitr::kable()
Month-Yr Dur priority_category_cd build_phase Day_of_Week vehicle_family_ref CPG_Uptime Month
Mar-2017:181 Durability :1122 P100:327 VP :522 Fri:160 JL :212 Min. :0.0000 Mar :181
Jul-2017:153 EBRG-Trans : 0 P70 :663 Proto :264 Mon:163 KL :130 1st Qu.:0.5288 Jul :153
Jun-2017:151 Reliability: 0 P85 :132 PS :135 Sat:152 M1 :103 Median :0.7770 Jun :151
May-2017:144 NA NA Mule :110 Sun:152 DT : 92 Mean :0.7023 May :144
Feb-2017:141 NA NA J1 : 30 Thu:163 MP : 83 3rd Qu.:0.9725 Feb :141
Apr-2017:124 NA NA Comp Veh: 21 Tue:166 RU-PHEV: 80 Max. :1.0000 Apr :124
(Other) :228 NA NA (Other) : 40 Wed:166 (Other):422 NA (Other):228 
CPGDurability=na.omit(CPGDurability)

CPGDurability_ts=ts(CPGDurability$CPG_Uptime)
# fitting weekly data. number of seasons is 52 

ss <- AddSemilocalLinearTrend(list(), Durability_ts)
ss <- AddSeasonal(ss, Durability_ts, nseasons = 52)
#model <- bsts(APG_Uptime~ ., state.specification = ss, niter = 1000,data =)
model <- bsts(Durability_ts, state.specification = ss, niter = 1000)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:05:57 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:06:01 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:06:06 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:06:10 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:06:15 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:06:19 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:06:24 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:06:29 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:06:34 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:06:38 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, model)
pred <- predict(model, horizon = 100, burn = burn, quantiles = c(.025, .975))
#pred
#pred
d2=data.frame(t(pred$interval),pred$mean,pred$median,Index=832:931)

names(d2) <- c("lower95", "upper95", "Actual","median","Index")

DT::datatable(d2) 
#model$one.step.prediction.errors #residual errors 831
#yhat=y-e
d3=data.frame(Actual=pred$original.series,Fitted=-colMeans(model$one.step.prediction.errors[-(1:burn),])+Durability_ts,Index=1:length(Durability_ts))


DT::datatable(d3)
#length(pred["bsts.prediction"])  #1



#apply(pred$distribution, 2, mean)



### Plot actual versus predicted with credible intervals for the holdout period
ggplot(data=d3, aes(x=Index)) +
   geom_line(aes(y=Actual, colour = "Actual"), size=1.5, linetype=1) +
   geom_line(aes(y=Fitted, colour = "Fitted"), size=1.5, linetype=1) +
   theme_bw() + theme(legend.title = element_blank()) +
   geom_line(aes(x = Index, y = Actual, colour="Forecast"), data =d2, alpha = 0.5, size=1.5, linetype=2) +
   geom_ribbon(aes(ymin=lower95, ymax=upper95), fill="grey", alpha=0.5,data=d2) +
   theme(axis.text.x=element_text(angle = -90, hjust = 0))+
  geom_smooth(aes(x = Index, y = Actual), data =d2,
                method = 'loess')+
geom_vline(xintercept=831, linetype=2) +
  scale_y_continuous(labels = percent_format())+
  labs(title = "Predicted Values vs Observed values vs Forecast", x = "Time Index",y="APG_Uptime")

#colMeans(pred$distribution) =pred$mean


d2=data.frame(Index=832:931,pred$mean,t(pred$interval),pred$median)
names(d2) <- c("Index", "Mean" , "lower95", "upper95","median")
rio::export(d2, "/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/Uptimebayes.xlsx")
plot(model)

plot(model, "components")  # plot(model1, "comp") works too!

plot(model, "help")
## starting httpd help server ... done
plot(pred, plot.original = 831) #plot predictions together with first 831 observations

Durability Vehicle Prediction

modeldata=read_excel("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/actvsfcdatajuly311.xlsx",sheet = "Sheet1")%>%rename(Date=Week)%>%na.omit()
modeldata%>%head()
dim(modeldata)
## [1] 86  2
modeldata%>%mutate(Date=as.Date(Date))%>%summary()
##       Date                Actual     
##  Min.   :2016-01-11   Min.   :21.00  
##  1st Qu.:2016-06-07   1st Qu.:27.00  
##  Median :2016-11-03   Median :30.00  
##  Mean   :2016-11-13   Mean   :29.57  
##  3rd Qu.:2017-04-22   3rd Qu.:32.00  
##  Max.   :2017-09-18   Max.   :40.00
modeldata%>%mutate(Date=as.Date(Date))%>%ggplot(aes(Date,Actual))+geom_line()+
  
  theme_tq()+ geom_point(color = palette_light()[[1]], alpha = 0.5) +
  labs(title = "Durability Vehicles Forecasting", x = "Time", y = "Number of \n Durability Vehicles",
       subtitle = "data from 2016 & 2017") +
  scale_y_continuous(labels = scales::comma) +
  theme_tq()+
  scale_x_date(labels = date_format("%m-%d-%Y"))

idx=modeldata%>%tk_index()%>%tk_get_timeseries_signature() 
idx
# Augmenting a data frame

md=modeldata%>%mutate(YQ = as.yearqtr(Date))
md=tk_augment_timeseries_signature(md)
md
# Convert tibble to ts object with tk_ts()

md1=md%>%dplyr::select(Date,Actual)
  
md_ts <- tk_ts(md1, silent = TRUE)

# get index from md tk_get_timeseries_signature()

md2=idx%>% tk_index()
md2%>%knitr::kable()
## Warning in kable_markdown(x = structure(c("2016-01-11", "2016-01-18",
## "2016-01-25", : The table should have a header (column names)
2016-01-11
2016-01-18
2016-01-25
2016-02-01
2016-02-08
2016-02-15
2016-02-22
2016-02-29
2016-03-07
2016-03-14
2016-03-21
2016-03-28
2016-04-04
2016-04-11
2016-04-18
2016-04-25
2016-05-02
2016-05-09
2016-05-16
2016-05-23
2016-05-30
2016-06-06
2016-06-13
2016-06-20
2016-06-27
2016-07-04
2016-07-11
2016-07-18
2016-07-25
2016-08-01
2016-08-08
2016-08-15
2016-08-22
2016-08-29
2016-09-05
2016-09-12
2016-09-19
2016-09-26
2016-10-03
2016-10-10
2016-10-17
2016-10-24
2016-10-31
2016-11-07
2016-11-14
2016-11-28
2016-12-05
2016-12-12
2017-01-02
2017-01-09
2017-01-16
2017-01-23
2017-01-30
2017-02-06
2017-02-13
2017-02-20
2017-02-27
2017-03-06
2017-03-13
2017-03-20
2017-03-27
2017-04-03
2017-04-10
2017-04-17
2017-04-24
2017-05-01
2017-05-08
2017-05-15
2017-05-22
2017-05-29
2017-06-05
2017-06-12
2017-06-19
2017-06-26
2017-07-03
2017-07-10
2017-07-17
2017-07-24
2017-07-31
2017-08-07
2017-08-14
2017-08-21
2017-08-28
2017-09-04
2017-09-11
2017-09-18 
idx_future<-md2 %>%
    tk_make_future_timeseries(n_future = 12)
idx_future
##  [1] "2017-09-25 UTC" "2017-10-02 UTC" "2017-10-09 UTC" "2017-10-16 UTC"
##  [5] "2017-10-23 UTC" "2017-10-30 UTC" "2017-11-06 UTC" "2017-11-13 UTC"
##  [9] "2017-11-20 UTC" "2017-11-27 UTC" "2017-12-04 UTC" "2017-12-11 UTC"
md%>%knitr::kable()
Date Actual YQ index.num diff year half quarter month month.xts month.lbl day hour minute second hour12 am.pm wday wday.xts wday.lbl mday qday yday mweek week week.iso week2 week3 week4 mday7
2016-01-11 26 2016 Q1 1452470400 NA 2016 1 1 1 0 January 11 0 0 0 0 1 2 1 Monday 11 11 11 3 2 2 0 2 2 2
2016-01-18 27 2016 Q1 1453075200 604800 2016 1 1 1 0 January 18 0 0 0 0 1 2 1 Monday 18 18 18 4 3 3 1 0 3 3
2016-01-25 28 2016 Q1 1453680000 604800 2016 1 1 1 0 January 25 0 0 0 0 1 2 1 Monday 25 25 25 5 4 4 0 1 0 4
2016-02-01 22 2016 Q1 1454284800 604800 2016 1 1 2 1 February 1 0 0 0 0 1 2 1 Monday 1 32 32 6 5 5 1 2 1 1
2016-02-08 27 2016 Q1 1454889600 604800 2016 1 1 2 1 February 8 0 0 0 0 1 2 1 Monday 8 39 39 2 6 6 0 0 2 2
2016-02-15 31 2016 Q1 1455494400 604800 2016 1 1 2 1 February 15 0 0 0 0 1 2 1 Monday 15 46 46 3 7 7 1 1 3 3
2016-02-22 26 2016 Q1 1456099200 604800 2016 1 1 2 1 February 22 0 0 0 0 1 2 1 Monday 22 53 53 4 8 8 0 2 0 4
2016-02-29 27 2016 Q1 1456704000 604800 2016 1 1 2 1 February 29 0 0 0 0 1 2 1 Monday 29 60 60 5 9 9 1 0 1 5
2016-03-07 21 2016 Q1 1457308800 604800 2016 1 1 3 2 March 7 0 0 0 0 1 2 1 Monday 7 67 67 2 10 10 0 1 2 2
2016-03-14 23 2016 Q1 1457913600 604800 2016 1 1 3 2 March 14 0 0 0 0 1 2 1 Monday 14 74 74 3 11 11 1 2 3 3
2016-03-21 25 2016 Q1 1458518400 604800 2016 1 1 3 2 March 21 0 0 0 0 1 2 1 Monday 21 81 81 4 12 12 0 0 0 4
2016-03-28 23 2016 Q1 1459123200 604800 2016 1 1 3 2 March 28 0 0 0 0 1 2 1 Monday 28 88 88 5 13 13 1 1 1 5
2016-04-04 23 2016 Q2 1459728000 604800 2016 1 2 4 3 April 4 0 0 0 0 1 2 1 Monday 4 4 95 2 14 14 0 2 2 1
2016-04-11 23 2016 Q2 1460332800 604800 2016 1 2 4 3 April 11 0 0 0 0 1 2 1 Monday 11 11 102 3 15 15 1 0 3 2
2016-04-18 23 2016 Q2 1460937600 604800 2016 1 2 4 3 April 18 0 0 0 0 1 2 1 Monday 18 18 109 4 16 16 0 1 0 3
2016-04-25 27 2016 Q2 1461542400 604800 2016 1 2 4 3 April 25 0 0 0 0 1 2 1 Monday 25 25 116 5 17 17 1 2 1 4
2016-05-02 31 2016 Q2 1462147200 604800 2016 1 2 5 4 May 2 0 0 0 0 1 2 1 Monday 2 32 123 1 18 18 0 0 2 1
2016-05-09 30 2016 Q2 1462752000 604800 2016 1 2 5 4 May 9 0 0 0 0 1 2 1 Monday 9 39 130 2 19 19 1 1 3 2
2016-05-16 33 2016 Q2 1463356800 604800 2016 1 2 5 4 May 16 0 0 0 0 1 2 1 Monday 16 46 137 3 20 20 0 2 0 3
2016-05-23 29 2016 Q2 1463961600 604800 2016 1 2 5 4 May 23 0 0 0 0 1 2 1 Monday 23 53 144 4 21 21 1 0 1 4
2016-05-30 29 2016 Q2 1464566400 604800 2016 1 2 5 4 May 30 0 0 0 0 1 2 1 Monday 30 60 151 5 22 22 0 1 2 5
2016-06-06 32 2016 Q2 1465171200 604800 2016 1 2 6 5 June 6 0 0 0 0 1 2 1 Monday 6 67 158 2 23 23 1 2 3 1
2016-06-13 31 2016 Q2 1465776000 604800 2016 1 2 6 5 June 13 0 0 0 0 1 2 1 Monday 13 74 165 3 24 24 0 0 0 2
2016-06-20 36 2016 Q2 1466380800 604800 2016 1 2 6 5 June 20 0 0 0 0 1 2 1 Monday 20 81 172 4 25 25 1 1 1 3
2016-06-27 37 2016 Q2 1466985600 604800 2016 1 2 6 5 June 27 0 0 0 0 1 2 1 Monday 27 88 179 5 26 26 0 2 2 4
2016-07-04 38 2016 Q3 1467590400 604800 2016 2 3 7 6 July 4 0 0 0 0 1 2 1 Monday 4 4 186 2 27 27 1 0 3 1
2016-07-11 39 2016 Q3 1468195200 604800 2016 2 3 7 6 July 11 0 0 0 0 1 2 1 Monday 11 11 193 3 28 28 0 1 0 2
2016-07-18 40 2016 Q3 1468800000 604800 2016 2 3 7 6 July 18 0 0 0 0 1 2 1 Monday 18 18 200 4 29 29 1 2 1 3
2016-07-25 35 2016 Q3 1469404800 604800 2016 2 3 7 6 July 25 0 0 0 0 1 2 1 Monday 25 25 207 5 30 30 0 0 2 4
2016-08-01 38 2016 Q3 1470009600 604800 2016 2 3 8 7 August 1 0 0 0 0 1 2 1 Monday 1 32 214 6 31 31 1 1 3 1
2016-08-08 32 2016 Q3 1470614400 604800 2016 2 3 8 7 August 8 0 0 0 0 1 2 1 Monday 8 39 221 2 32 32 0 2 0 2
2016-08-15 33 2016 Q3 1471219200 604800 2016 2 3 8 7 August 15 0 0 0 0 1 2 1 Monday 15 46 228 3 33 33 1 0 1 3
2016-08-22 29 2016 Q3 1471824000 604800 2016 2 3 8 7 August 22 0 0 0 0 1 2 1 Monday 22 53 235 4 34 34 0 1 2 4
2016-08-29 29 2016 Q3 1472428800 604800 2016 2 3 8 7 August 29 0 0 0 0 1 2 1 Monday 29 60 242 5 35 35 1 2 3 5
2016-09-05 27 2016 Q3 1473033600 604800 2016 2 3 9 8 September 5 0 0 0 0 1 2 1 Monday 5 67 249 2 36 36 0 0 0 1
2016-09-12 24 2016 Q3 1473638400 604800 2016 2 3 9 8 September 12 0 0 0 0 1 2 1 Monday 12 74 256 3 37 37 1 1 1 2
2016-09-19 23 2016 Q3 1474243200 604800 2016 2 3 9 8 September 19 0 0 0 0 1 2 1 Monday 19 81 263 4 38 38 0 2 2 3
2016-09-26 31 2016 Q3 1474848000 604800 2016 2 3 9 8 September 26 0 0 0 0 1 2 1 Monday 26 88 270 5 39 39 1 0 3 4
2016-10-03 32 2016 Q4 1475452800 604800 2016 2 4 10 9 October 3 0 0 0 0 1 2 1 Monday 3 3 277 2 40 40 0 1 0 1
2016-10-10 33 2016 Q4 1476057600 604800 2016 2 4 10 9 October 10 0 0 0 0 1 2 1 Monday 10 10 284 3 41 41 1 2 1 2
2016-10-17 30 2016 Q4 1476662400 604800 2016 2 4 10 9 October 17 0 0 0 0 1 2 1 Monday 17 17 291 4 42 42 0 0 2 3
2016-10-24 32 2016 Q4 1477267200 604800 2016 2 4 10 9 October 24 0 0 0 0 1 2 1 Monday 24 24 298 5 43 43 1 1 3 4
2016-10-31 33 2016 Q4 1477872000 604800 2016 2 4 10 9 October 31 0 0 0 0 1 2 1 Monday 31 31 305 6 44 44 0 2 0 5
2016-11-07 35 2016 Q4 1478476800 604800 2016 2 4 11 10 November 7 0 0 0 0 1 2 1 Monday 7 38 312 2 45 45 1 0 1 2
2016-11-14 32 2016 Q4 1479081600 604800 2016 2 4 11 10 November 14 0 0 0 0 1 2 1 Monday 14 45 319 3 46 46 0 1 2 3
2016-11-28 30 2016 Q4 1480291200 1209600 2016 2 4 11 10 November 28 0 0 0 0 1 2 1 Monday 28 59 333 5 48 48 0 0 0 5
2016-12-05 29 2016 Q4 1480896000 604800 2016 2 4 12 11 December 5 0 0 0 0 1 2 1 Monday 5 66 340 2 49 49 1 1 1 1
2016-12-12 27 2016 Q4 1481500800 604800 2016 2 4 12 11 December 12 0 0 0 0 1 2 1 Monday 12 73 347 3 50 50 0 2 2 2
2017-01-02 27 2017 Q1 1483315200 1814400 2017 1 1 1 0 January 2 0 0 0 0 1 2 1 Monday 2 2 2 1 1 1 1 1 1 1
2017-01-09 26 2017 Q1 1483920000 604800 2017 1 1 1 0 January 9 0 0 0 0 1 2 1 Monday 9 9 9 2 2 2 0 2 2 2
2017-01-16 26 2017 Q1 1484524800 604800 2017 1 1 1 0 January 16 0 0 0 0 1 2 1 Monday 16 16 16 3 3 3 1 0 3 3
2017-01-23 25 2017 Q1 1485129600 604800 2017 1 1 1 0 January 23 0 0 0 0 1 2 1 Monday 23 23 23 4 4 4 0 1 0 4
2017-01-30 25 2017 Q1 1485734400 604800 2017 1 1 1 0 January 30 0 0 0 0 1 2 1 Monday 30 30 30 5 5 5 1 2 1 5
2017-02-06 22 2017 Q1 1486339200 604800 2017 1 1 2 1 February 6 0 0 0 0 1 2 1 Monday 6 37 37 2 6 6 0 0 2 1
2017-02-13 25 2017 Q1 1486944000 604800 2017 1 1 2 1 February 13 0 0 0 0 1 2 1 Monday 13 44 44 3 7 7 1 1 3 2
2017-02-20 26 2017 Q1 1487548800 604800 2017 1 1 2 1 February 20 0 0 0 0 1 2 1 Monday 20 51 51 4 8 8 0 2 0 3
2017-02-27 23 2017 Q1 1488153600 604800 2017 1 1 2 1 February 27 0 0 0 0 1 2 1 Monday 27 58 58 5 9 9 1 0 1 4
2017-03-06 25 2017 Q1 1488758400 604800 2017 1 1 3 2 March 6 0 0 0 0 1 2 1 Monday 6 65 65 2 10 10 0 1 2 1
2017-03-13 28 2017 Q1 1489363200 604800 2017 1 1 3 2 March 13 0 0 0 0 1 2 1 Monday 13 72 72 3 11 11 1 2 3 2
2017-03-20 28 2017 Q1 1489968000 604800 2017 1 1 3 2 March 20 0 0 0 0 1 2 1 Monday 20 79 79 4 12 12 0 0 0 3
2017-03-27 29 2017 Q1 1490572800 604800 2017 1 1 3 2 March 27 0 0 0 0 1 2 1 Monday 27 86 86 5 13 13 1 1 1 4
2017-04-03 30 2017 Q2 1491177600 604800 2017 1 2 4 3 April 3 0 0 0 0 1 2 1 Monday 3 3 93 2 14 14 0 2 2 1
2017-04-10 32 2017 Q2 1491782400 604800 2017 1 2 4 3 April 10 0 0 0 0 1 2 1 Monday 10 10 100 3 15 15 1 0 3 2
2017-04-17 36 2017 Q2 1492387200 604800 2017 1 2 4 3 April 17 0 0 0 0 1 2 1 Monday 17 17 107 4 16 16 0 1 0 3
2017-04-24 32 2017 Q2 1492992000 604800 2017 1 2 4 3 April 24 0 0 0 0 1 2 1 Monday 24 24 114 5 17 17 1 2 1 4
2017-05-01 27 2017 Q2 1493596800 604800 2017 1 2 5 4 May 1 0 0 0 0 1 2 1 Monday 1 31 121 6 18 18 0 0 2 1
2017-05-08 32 2017 Q2 1494201600 604800 2017 1 2 5 4 May 8 0 0 0 0 1 2 1 Monday 8 38 128 2 19 19 1 1 3 2
2017-05-15 28 2017 Q2 1494806400 604800 2017 1 2 5 4 May 15 0 0 0 0 1 2 1 Monday 15 45 135 3 20 20 0 2 0 3
2017-05-22 28 2017 Q2 1495411200 604800 2017 1 2 5 4 May 22 0 0 0 0 1 2 1 Monday 22 52 142 4 21 21 1 0 1 4
2017-05-29 30 2017 Q2 1496016000 604800 2017 1 2 5 4 May 29 0 0 0 0 1 2 1 Monday 29 59 149 5 22 22 0 1 2 5
2017-06-05 34 2017 Q2 1496620800 604800 2017 1 2 6 5 June 5 0 0 0 0 1 2 1 Monday 5 66 156 2 23 23 1 2 3 1
2017-06-12 32 2017 Q2 1497225600 604800 2017 1 2 6 5 June 12 0 0 0 0 1 2 1 Monday 12 73 163 3 24 24 0 0 0 2
2017-06-19 28 2017 Q2 1497830400 604800 2017 1 2 6 5 June 19 0 0 0 0 1 2 1 Monday 19 80 170 4 25 25 1 1 1 3
2017-06-26 34 2017 Q2 1498435200 604800 2017 1 2 6 5 June 26 0 0 0 0 1 2 1 Monday 26 87 177 5 26 26 0 2 2 4
2017-07-03 32 2017 Q3 1499040000 604800 2017 2 3 7 6 July 3 0 0 0 0 1 2 1 Monday 3 3 184 2 27 27 1 0 3 1
2017-07-10 35 2017 Q3 1499644800 604800 2017 2 3 7 6 July 10 0 0 0 0 1 2 1 Monday 10 10 191 3 28 28 0 1 0 2
2017-07-17 32 2017 Q3 1500249600 604800 2017 2 3 7 6 July 17 0 0 0 0 1 2 1 Monday 17 17 198 4 29 29 1 2 1 3
2017-07-24 30 2017 Q3 1500854400 604800 2017 2 3 7 6 July 24 0 0 0 0 1 2 1 Monday 24 24 205 5 30 30 0 0 2 4
2017-07-31 33 2017 Q3 1501459200 604800 2017 2 3 7 6 July 31 0 0 0 0 1 2 1 Monday 31 31 212 6 31 31 1 1 3 5
2017-08-07 33 2017 Q3 1502064000 604800 2017 2 3 8 7 August 7 0 0 0 0 1 2 1 Monday 7 38 219 2 32 32 0 2 0 2
2017-08-14 32 2017 Q3 1502668800 604800 2017 2 3 8 7 August 14 0 0 0 0 1 2 1 Monday 14 45 226 3 33 33 1 0 1 3
2017-08-21 30 2017 Q3 1503273600 604800 2017 2 3 8 7 August 21 0 0 0 0 1 2 1 Monday 21 52 233 4 34 34 0 1 2 4
2017-08-28 27 2017 Q3 1503878400 604800 2017 2 3 8 7 August 28 0 0 0 0 1 2 1 Monday 28 59 240 5 35 35 1 2 3 5
2017-09-04 32 2017 Q3 1504483200 604800 2017 2 3 9 8 September 4 0 0 0 0 1 2 1 Monday 4 66 247 2 36 36 0 0 0 1
2017-09-11 32 2017 Q3 1505088000 604800 2017 2 3 9 8 September 11 0 0 0 0 1 2 1 Monday 11 73 254 3 37 37 1 1 1 2
2017-09-18 36 2017 Q3 1505692800 604800 2017 2 3 9 8 September 18 0 0 0 0 1 2 1 Monday 18 80 261 4 38 38 0 2 2 3 
# fitting weekly data. number of seasons is 52 

ss <- AddSemilocalLinearTrend(list(), md_ts)
ss <- AddSeasonal(ss, md_ts, nseasons = 52)
modeld <- bsts(md_ts, state.specification = ss, niter = 1000)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:07:57 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:07:57 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:07:58 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:07:58 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:07:59 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:07:59 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:08:00 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:08:00 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:08:01 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:08:01 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, modeld)
predd <- predict(modeld, horizon = 12, burn = burn, quantiles = c(.025, .975))
#predd$mean

d4=data.frame(idx_future,predd$mean,t(predd$interval),predd$median)
names(d4) <- c("Time", "Mean" , "lower95", "upper95","median")
d4%>%knitr::kable()
Time Mean lower95 upper95 median
2017-09-25 33.50818 25.98688 41.31687 33.52783
2017-10-02 31.71297 22.18931 40.64716 31.85725
2017-10-09 30.25796 20.17723 40.74714 30.15724
2017-10-16 36.52221 25.06697 47.55274 36.27195
2017-10-23 36.52045 23.38043 49.15911 36.34008
2017-10-30 36.90029 23.81006 50.05009 36.84632
2017-11-06 34.19759 19.71037 48.28863 34.41766
2017-11-13 35.59264 19.83561 50.17114 35.59512
2017-11-20 36.43044 19.04633 51.37173 36.96453
2017-11-27 38.13852 19.70805 53.85161 38.49782
2017-12-04 36.43115 17.64673 53.74799 36.59691
2017-12-11 35.42211 15.94861 53.68617 35.65926 
rio::export(d4, "/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/Durabilitybayes.xlsx")
#model$one.step.prediction.errors #residual errors 831
#yhat=y-e
d5=data.frame(Time=modeldata$Date,Actual=predd$original.series,Fit=-colMeans(modeld$one.step.prediction.errors[-(1:burn),])+md_ts)%>%rename(Fit=Actual.1)


DT::datatable(d5)
#length(pred["bsts.prediction"])  #1



#apply(pred$distribution, 2, mean)

names(d4) <- c("Time", "Actual" , "lower95", "upper95","median")

DT::datatable(d4)
### Plot actual versus predicted with credible intervals for the holdout period
ggplot(data=d5, aes(x=as.Date(Time))) +
   geom_line(aes(y=Actual, colour = "Actual Data"), size=1.5, linetype=1) +
   geom_line(aes(y=Fit, colour = "Model Fit on Data"), size=1.5, linetype=1) +
   theme_bw() + theme(legend.title = element_blank()) +
  
   geom_ribbon(aes(ymin=lower95, ymax=upper95), fill="grey", alpha=0.5,data=d4) +

  theme(axis.text.x=element_text(angle=50, size=10, vjust=0.5))+
 
 geom_vline(xintercept=as.Date("2017-09-18"), linetype=2)+

  labs(title = "Predicted Values vs Observed values vs Forecast\n  Bayesian Structural Model", x = "Time ",y="Durability Vehicles")+
geom_line(aes(x =as.Date(Time), y = Actual, colour="12-Week Forecast"), data =d4, alpha = 0.5, size=1.5, linetype=1) +


(scale_x_date(breaks=date_breaks("1 month"),
      labels=date_format("%b %d  %y")))

### Plot actual versus predicted with credible intervals for the holdout period
ggplot(data=d5, aes(x=as.Date(Time))) +
   geom_line(aes(y=Actual, colour = "Actual"), size=1.5, linetype=1) +
   geom_line(aes(y=Fit, colour = "Fit"), size=1.5, linetype=1) +
   theme_bw() + theme(legend.title = element_blank()) +
  
   geom_ribbon(aes(ymin=lower95, ymax=upper95), fill="grey", alpha=0.5,data=d4) +
   #theme(axis.text.x=element_text(angle = 90, hjust = 0))+
  theme(axis.text.x=element_text(angle=50, size=10, vjust=0.5))+
 
 geom_vline(xintercept=as.Date("2017-09-18"), linetype=2)+

  labs(title = "Predicted Values vs Observed values vs Forecast\n  Bayesian Structural Model", x = "Time ",y="Durability Vehicles")+
geom_line(aes(x =as.Date(Time), y = Actual, colour="Forecast"), data =d4, alpha = 0.5, size=1.5, linetype=1) +
 geom_smooth(aes(x = as.Date(Time), y = Actual), data =d4,
                method = 'loess')+

(scale_x_date(breaks=date_breaks("1 month"),
      labels=date_format("%b %d  %y")))

xreg=Durability%>%select(Month,vehicle_family_ref,Day_of_Week,Build_Category,build_phase,priority_category_cd)

xreg2=model.matrix(~Month+vehicle_family_ref+Day_of_Week+Build_Category+build_phase+priority_category_cd, xreg)

xreg3=Matrix::sparse.model.matrix(~Month+vehicle_family_ref+Day_of_Week+Build_Category+build_phase, data.frame(xreg))

xreg4=Matrix::sparse.model.matrix(~Month+vehicle_family_ref+Build_Category, data.frame(xreg))



xreg5=model.matrix(~Month+vehicle_family_ref, data.frame(xreg))

head(xreg4)
## [1] 1 1 1 1 1 1
# fitting weekly data. number of seasons is 52 

ss <- AddSemilocalLinearTrend(list(), Durability_ts)
ss <- AddSeasonal(ss,Durability_ts, nseasons = 52)


model4 <- bsts(APG_Uptime~ Month+vehicle_family_ref+Build_Category+Day_of_Week, state.specification = ss, niter = 1000,data =Durability)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:08:05 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:08:10 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:08:15 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:08:20 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:08:25 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:08:30 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:08:34 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:08:39 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:08:44 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:08:49 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, modeld)




#preda <- predict(model4, horizon = 100, burn = burn, quantiles = c(.025, .975))
#Error in .ExtractPredictors(object, newdata, na.action = na.action) : 
#  You need to supply 'newdata' when making predictions with a bsts object that has a #regression component.
model4$coefficients%>%head()%>%knitr::kable()
(Intercept) MonthAug MonthFeb MonthJan MonthJul MonthJun MonthMar MonthMay vehicle_family_refD2 vehicle_family_refDS vehicle_family_refDT vehicle_family_refJK vehicle_family_refJL vehicle_family_refJT vehicle_family_refKL vehicle_family_refLA-SRT vehicle_family_refLD vehicle_family_refLX-SRT vehicle_family_refMP vehicle_family_refRU vehicle_family_refRU-PHEV vehicle_family_refVM vehicle_family_refWD-SRT vehicle_family_refWK vehicle_family_refWK-SRT Build_CategoryJ1 Build_CategoryMule_Proto Build_CategoryPS Build_CategoryVP Day_of_WeekMon Day_of_WeekSat Day_of_WeekSun Day_of_WeekThu Day_of_WeekTue Day_of_WeekWed
0 0 0 0 0 0 0 0 0 0 0 0 -0.1126601 0 0 0 0 0 0 0 0 0 0 0 -0.3339215 0.4003750 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 -0.1089064 0 0 0 0 0 0 0 0 0 0 0 -0.2171084 0.4349716 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 -0.1269382 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.4534998 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.5060575 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.5598902 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.6137876 0 0 0 0 0 0 0 0 0 
model4$state.contributions%>%head()%>%knitr::kable()
## Warning in kable_markdown(x = structure(c("0.0338861", "0.0481932",
## "-0.0535962", : The table should have a header (column names)
0.0338861
0.0481932
-0.0535962
-0.1016431
0.1009437
0.2550663 
model4$model.options%>%head()%>%knitr::kable()
## Warning in kable_markdown(x = structure("TRUE", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
## Warning in kable_markdown(x = structure("TRUE", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
## Warning in kable_markdown(x = structure("SSVS", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
## Warning in kable_markdown(x = structure("0", .Dim = c(1L, 1L), .Dimnames =
## list(: The table should have a header (column names)
## Warning in kable_markdown(x = structure("0.01", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
## Warning in kable_markdown(x = structure("Inf", .Dim = c(1L, 1L), .Dimnames
## = list(: The table should have a header (column names)
TRUE
TRUE
SSVS
0
0.01
Inf 
model4$predictors%>%head()%>%knitr::kable()
(Intercept) MonthAug MonthFeb MonthJan MonthJul MonthJun MonthMar MonthMay vehicle_family_refD2 vehicle_family_refDS vehicle_family_refDT vehicle_family_refJK vehicle_family_refJL vehicle_family_refJT vehicle_family_refKL vehicle_family_refLA-SRT vehicle_family_refLD vehicle_family_refLX-SRT vehicle_family_refMP vehicle_family_refRU vehicle_family_refRU-PHEV vehicle_family_refVM vehicle_family_refWD-SRT vehicle_family_refWK vehicle_family_refWK-SRT Build_CategoryJ1 Build_CategoryMule_Proto Build_CategoryPS Build_CategoryVP Day_of_WeekMon Day_of_WeekSat Day_of_WeekSun Day_of_WeekThu Day_of_WeekTue Day_of_WeekWed
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
model4$terms%>%head()
## APG_Uptime ~ Month + vehicle_family_ref + Build_Category
## attr(,"variables")
## list(APG_Uptime, Month, vehicle_family_ref, Build_Category)
## attr(,"factors")
##                    Month vehicle_family_ref Build_Category
## APG_Uptime             0                  0              0
## Month                  1                  0              0
## vehicle_family_ref     0                  1              0
## Build_Category         0                  0              1
## attr(,"term.labels")
## [1] "Month"              "vehicle_family_ref" "Build_Category"    
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(APG_Uptime, Month, vehicle_family_ref, Build_Category)
## attr(,"dataClasses")
##         APG_Uptime              Month vehicle_family_ref 
##          "numeric"           "factor"           "factor" 
##     Build_Category 
##           "factor"
model4$state.contributions%>%head()%>%knitr::kable()
## Warning in kable_markdown(x = structure(c("0.0338861", "0.0481932",
## "-0.0535962", : The table should have a header (column names)
0.0338861
0.0481932
-0.0535962
-0.1016431
0.1009437
0.2550663 
plot(model4, "coef")

plot(model4, "components")  # plot(model1, "comp") works too!

ss <- AddSemilocalLinearTrend(list(), Durability_ts)
ss <- AddSeasonal(ss,Durability_ts, nseasons = 52)


model6 <- bsts(APG_Uptime~ vehicle_family_ref+Build_Category, state.specification = ss, niter = 1000,data =Durability)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:10:11 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:10:16 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:10:21 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:10:25 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:10:30 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:10:35 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:10:39 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:10:44 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:10:49 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:10:53 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, modeld)

plot(model6, "coef")

ss <- AddSemilocalLinearTrend(list(), Durability_ts)
ss <- AddSeasonal(ss,Durability_ts, nseasons = 52)


model5 <- bsts(Durability_ts, state.specification = ss, niter = 1000,data =Durability)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:10:58 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:11:03 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:11:07 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:11:12 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:11:16 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:11:21 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:11:25 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:11:29 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:11:34 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:11:38 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, modeld)

m4=bsts.prediction.errors(model4)

m5=bsts.prediction.errors(model5)


CompareBstsModels(list("Model 1" = model5,
                       "Model 2" = model4),
                  colors = c("red", "blue"))

The addition of the regressors did not improve the model prediction as there is no difference between the model with regressors and the response variable and the model with only the response variable.

# fitting weekly data. number of seasons is 52 

ss <- AddSemilocalLinearTrend(list(), CPGDurability_ts)
ss <- AddSeasonal(ss, CPGDurability_ts, nseasons = 52)
#model <- bsts(APG_Uptime~ ., state.specification = ss, niter = 1000,data =)
modelCPG <- bsts(CPGDurability_ts, state.specification = ss, niter = 1000)
## =-=-=-=-= Iteration 0 Mon Sep 18 22:11:44 2017 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Sep 18 22:11:50 2017 =-=-=-=-=
## =-=-=-=-= Iteration 200 Mon Sep 18 22:11:56 2017 =-=-=-=-=
## =-=-=-=-= Iteration 300 Mon Sep 18 22:12:02 2017 =-=-=-=-=
## =-=-=-=-= Iteration 400 Mon Sep 18 22:12:08 2017 =-=-=-=-=
## =-=-=-=-= Iteration 500 Mon Sep 18 22:12:14 2017 =-=-=-=-=
## =-=-=-=-= Iteration 600 Mon Sep 18 22:12:20 2017 =-=-=-=-=
## =-=-=-=-= Iteration 700 Mon Sep 18 22:12:27 2017 =-=-=-=-=
## =-=-=-=-= Iteration 800 Mon Sep 18 22:12:33 2017 =-=-=-=-=
## =-=-=-=-= Iteration 900 Mon Sep 18 22:12:39 2017 =-=-=-=-=
burn <- SuggestBurn(0.1, model)
predCPG <- predict(modelCPG, horizon = 100, burn = burn, quantiles = c(.025, .975))
plot(predCPG, plot.original = 1122) #plot predictions together with first 831 observations

plot(modelCPG)

plot(model, "components")  # plot(model1, "comp") works too!

d2CPG=data.frame(t(predCPG$interval),predCPG$mean,predCPG$median,Index=1123:1222)

names(d2CPG) <- c("lower95", "upper95", "Actual","median","Index")

DT::datatable(d2CPG) 
#model$one.step.prediction.errors #residual errors 831
#yhat=y-e
d3CPG=data.frame(Actual=predCPG$original.series,Fitted=-colMeans(modelCPG$one.step.prediction.errors[-(1:burn),])+CPGDurability_ts,Index=1:length(CPGDurability_ts))


DT::datatable(d3CPG)
#length(pred["bsts.prediction"])  #1



#apply(pred$distribution, 2, mean)



### Plot actual versus predicted with credible intervals for the holdout period
ggplot(data=d3CPG, aes(x=Index)) +
   geom_line(aes(y=Actual, colour = "Actual"), size=1.5, linetype=1) +
   geom_line(aes(y=Fitted, colour = "Fitted"), size=1.5, linetype=1) +
   theme_bw() + theme(legend.title = element_blank()) +
   geom_point(aes(x = Index, y = Actual, colour="Forecast"), data =d2CPG, alpha = 0.5, size=1.5, linetype=2) +
   geom_ribbon(aes(ymin=lower95, ymax=upper95), fill="grey", alpha=0.5,data=d2CPG) +
   theme(axis.text.x=element_text(angle = -90, hjust = 0))+
  #geom_smooth(aes(x = Index, y = Actual), data =d2CPG,method = 'loess')+
geom_vline(xintercept=1122, linetype=2) +
  scale_y_continuous(labels = percent_format())+
  labs(title = "Predicted Values vs Observed values vs Forecast\n Bayesian Structural Model", x = "Time Index",y="CPG Uptime")
## Warning: Ignoring unknown parameters: linetype

d2CPG=data.frame(Index=1123:1222,t(predCPG$interval),predCPG$mean,predCPG$median)

names(d2CPG) <- c("Index","lower95", "upper95", "Actual","median")

DT::datatable(d2CPG) 
#model$one.step.prediction.errors #residual errors 831
#yhat=y-e
d3CPG=data.frame(Actual=predCPG$original.series,Fitted=-colMeans(modelCPG$one.step.prediction.errors[-(1:burn),])+CPGDurability_ts,Index=1:length(CPGDurability_ts))


DT::datatable(d3CPG)
#length(pred["bsts.prediction"])  #1



#apply(pred$distribution, 2, mean)



### Plot actual versus predicted with credible intervals for the holdout period
ggplot(data=d3CPG, aes(x=Index)) +
   geom_line(aes(y=Actual, colour = "Actual"), size=1.5, linetype=1) +
   geom_line(aes(y=Fitted, colour = "Fitted"), size=1.5, linetype=1) +
   theme_bw() + theme(legend.title = element_blank()) +
   geom_point(aes(x = Index, y = Actual, colour="Forecast"), data =d2CPG, alpha = 0.5, size=1.5, linetype=2) +
   geom_ribbon(aes(ymin=lower95, ymax=upper95), fill="grey", alpha=0.5,data=d2CPG) +
   theme(axis.text.x=element_text(angle = -90, hjust = 0))+
  geom_smooth(aes(x = Index, y = Actual), data =d2CPG,method = 'loess')+
geom_vline(xintercept=1122, linetype=2) +
  scale_y_continuous(labels = percent_format())+
  labs(title = "Predicted Values vs Observed values vs Forecast\n Bayesian Structural Model", x = "Time Index",y="CPG Uptime")
## Warning: Ignoring unknown parameters: linetype

names(d2CPG) <- c("Index","lower95", "upper95", "Forecast","median")
rio::export(d2CPG, "/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/Time series/Manpower/CPGDurabilitybayes.xlsx")