#CHAPTER 10: THROUGH THE LENS OF TIME SERIES
####################################################

####################################################
#LIBRARIES TO LOAD
####################################################
library(xlsx); library(ggfortify); library(forecast)
library(lubridate)

#Instructions for ggplot2-
#https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_ts.html

####################################################
##LOAD DATA
####################################################
#Import tblMainEmploy and name MainMon
#<path> provided by the user to inform of location of data
MaineMon<- read.xlsx(
	"C:\\<path>\\SkillsForCareerSecurity_Datasets.xlsx", 
	sheetName="tblMaineEmploy", header=TRUE) 

####################################################
####################################################
##10.1. WORKING WITH TIME SERIES
####################################################
##10.1.1. TIME SERIES AS AN R OBJECT
####################################################
	#BEGIN TS OBJECT

#AirPassengers dataset from R assigned to AP
data(AirPassengers)
AP<- AirPassengers
AP

#Characteristics of the object
class(AP)
start(AP)
end(AP)
frequency(AP)

#Charted series, year aggregated as trend, cycles
layout(1:3)
plot(AP, ylab="Passengers (1000's)", lwd=2)
plot(aggregate(AP), lwd=2)
boxplot(AP~cycle(AP), lwd=2)

	#END TS OBJECT
####################################################
####################################################
##10.1.2. TRANSFORMING DATA TO A TS OBJECT
####################################################
	#BEGIN TRANSFORMING TO TS OBJECT

#Import tblMainEmploy and name MainMon
MaineMon<- read.xlsx(navPath, 
	sheetName="tblMaineEmploy", header=TRUE) 

#Format MainMon to time series object
MaineMonTs<- ts(MaineMon$unemploy, start=c(1996,1),freq=12)
MaineMonTs

#Create aggregate series as annual
MaineAnnTs<- ts(aggregate(MaineMonTs, FUN = mean))
MaineAnnTs

#Chart monthly and annual contrast
layout(1:2) 
plot(MaineMonTs, ylab="Maine unemployed (%)")
plot(MaineAnnTs, ylab="Maine unemployed (%)")

	#END TRANSFORMING TO TS OBJECT	
####################################################
##10.1.3. WINDOWS TO TIME SERIES
####################################################
	#BEGIN WINDOWS

#Computation with time series
Maine.Feb<- window(MaineMonTs, start=c(1996,2), freq=T)  
Maine.Aug<- window(MaineMonTs, start=c(1996,8), freq=T)
(Feb.ratio<- mean(Maine.Feb)/mean(MaineMonTs) )
(Aug.ratio<- mean(Maine.Aug)/mean(MaineMonTs) )

#Import US unemployment data
UsEmploy<- read.xlsx(navPath, 
	sheetName="tblUsEmploy", header=TRUE) 

#Format as a time series
UsMonTs<- ts(UsEmploy$USun, start=c(1996,1), freq=12)

#Comparison of Maine to US unemployment
layout(1:3) 
plot(MaineMonTs, ylab="Maine unemployed (%)", lwd=2)
plot(MaineAnnTs, ylab="Maine unemployed (%)", lwd=2)
plot(UsMonTs, ylab="US unemployed(%)", lwd=2)

	#END WINDOWS
####################################################
##10.1.4. INTERSECTION OF SERIES
####################################################
	#BEGIN INTERSECTION

#Import a dataset of sequential data
#<path> provided by the user to inform of location of data
ChBrEl<- read.xlsx(
	<path>\\SkillsForCareerSecurity_Datasets.xlsx",
	sheetName="tblChocBeerElect", header=TRUE) 
head(ChBrEl)

#Convert the variables to the dataset to ts objects
Elec.ts<- ts(ChBrEl[,3], start=1958, freq=12)
Beer.ts<- ts(ChBrEl[,2], start=1958, freq=12)
Choc.ts<- ts(ChBrEl[,1], start=1958, freq=12)
plot(cbind(Elec.ts,Beer.ts,Choc.ts), lwd=2)
start(Elec.ts); end(Elec.ts)

#Import AirPassengers ts
data(AirPassengers)
AP<- AirPassengers 

#Create ts object for the mutual periods of data
ApElec<- ts.intersect(AP, Elec.ts)

start(ApElec)
end(ApElec)

#Plot of ApElec object.
layout(1:1)
plot(ApElec, main="Passengers and ELectrical", lwd=2)

	#END INTERSECTION
####################################################
####################################################
##10.2. ANALYTICS WITH TIME SERIES
####################################################
##10.2.2. DECOMPOSITION OF THE SERIES
####################################################	
	#BEGIN DECOMPOSITION

#For additive case
Elec.decom.add<- decompose(Elec.ts, type="add")
plot(Elec.decom.add, lwd=2)  

Elec.decom.mult<- decompose(Elec.ts, type="mult")
plot(Elec.decom.mult, lwd=2)

str(Elec.decom.add)
str(Elec.decom.mult)

#Trend distinguished from cycle
layout(1:1)
Trend<- Elec.decom.mult$trend
Seasonal<- Elec.decom.mult$seasonal
ts.plot(cbind(Trend, Trend*Seasonal), lty=1:2, lwd=2)

	#END DECOMPOSITION
####################################################
##10.2.3. AUTOCORRELATION AND CORRELOGRAMS
####################################################
	#BEGIN AUTOCORRELATION/CORRELOGRAMS

#Load dataset AirPassengers
data(AirPassengers)
AP<- AirPassengers

#Correlogram of observations
acf(AP, lwd=2)	

#Decompose AP and get count of periods
AP.decom<- decompose(AP, "multiplicative")
length(AP.decom$trend)

#Plot the random component and correlogram
layout(1:2)
plot(ts(AP.decom$random[7:138]), lwd=2)
acf(AP.decom$random[7:138], lwd=2)

#Correlogram with ARIMA
d.arima <- auto.arima(AirPassengers)
acf(d.arima$residuals, lwd=2)

	#END AUTOCORRELATION/CORRELOGRAMS
####################################################
##10.2.4. LEAD-LAG VARIABLES
####################################################
	#BEGIN LEAD-LAG

#Load data
Build<- read.xlsx(navPath,
	sheetName="tblApprvBuild", header=TRUE) 
head(Build)

#Plot ts objects in same chart
layout(1:1)
App.ts<- ts(Build$Approvals, start=c(1996,1), freq=4)
Act.ts<- ts(Build$Activity, start=c(1996,1), freq=4)
ts.plot(App.ts, Act.ts, lty=c(1,3), lwd=2)
abline(v=2000.5, lwd=3)

##Extract random components for approve and activity
#Approvals
app.ran<- decompose(App.ts)$random
app.ran.ts<- window(app.ran, start=c(1996,3), end=c(2006,1))
#Activity
act.ran<- decompose(Act.ts)$random
act.ran.ts<- window(act.ran, start=c(1996,3), end=c(2006,1)) ##Start period 3 due to NA

#Plot cross correlograms for observations and random
#Cross correlation of observations
layout(1:2)
ccf(App.ts, Act.ts, main="Observation Activity Lag Approvals")
#Cross correlogram of random component
ccf(app.ran.ts, act.ran.ts, main="Random Activity Lag Approvals")

	#END LEAD-LAG
####################################################
##10.2.5. SEVEN-DAY CYCLES
####################################################
	#BEGIN SEVEN-DAY

#Import dataset
#<path> provided by the user to inform of location of data
SevenDay<- read.xlsx(
	"C:\\<path>\\SkillsForCareerSecurity_Datasets.xlsx"",
	sheetName="tblWkCycle", header=TRUE) 
head(SevenDay)

#Transform series variable to ts object
WeekTs<- ts(SevenDay$Series, start=c(1,1), freq=7)

#Pattern for week
layout(1:3)
plot(WeekTs, lwd=2)
plot(aggregate(WeekTs), lwd=2)
boxplot(WeekTs~cycle(WeekTs))

#Decompose weeks as trend: Additive
DecomPredAdd<- decompose(WeekTs, type="add")
plot(DecomPredAdd, lwd=2)

#Decompose weeks as trend: Multiplicative
DecomPredMult<- decompose(WeekTs, type="mult")
plot(DecomPredMult, lwd=2)

#Extract seasonal component from respective decompositions
WkCycleAdd<- as.numeric(window(DecomPredAdd$seasonal, start=1, end=(1+6/7)))
WkCycleMult<- as.numeric(window(DecomPredMult$seasonal, start=1, end=(1+6/7)))

#Form data frame of additive and multiplicative
WeekDf<- data.frame(WkCycleAdd, WkCycleMult, DayNum = c(1:7), 
	DayTitle = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))
WeekDf

	#END SEVEN-DAY
####################################################
####################################################
##10.3. FORECASTING THE FUTURE
####################################################
##10.3.1. HOLT-WINTERS
####################################################
	#BEGIN HOLT-WINTERS METHOD

#Import AirPassengers dataset
data(AirPassengers)
AP<- AirPassengers

#Model with HoltWinters
AP.hw<- HoltWinters(AP, seasonal="mult")
AP.hw

#Plots for insight
plot(AP.hw, lwd=2)
plot(AP.hw$fitted, lwd=2)

	#END HOLT-WINTERS METHOD
####################################################
	#BEGIN FORECAST WITH HOLT-WINTERS

#Forecast with Holt Winters
#Predict function to forecast 4 years ahead
AP.predict<- predict(AP.hw, n.ahead=4*12)
#Plot forecast from past
ts.plot(AP, AP.predict, lty=1:2, col=c("black", "red"), lwd=2:3)

	#END FORECAST HOLT-WINTERS
####################################################
	#BEGIN TABLE OF FORECAST AND PLOT

#Convert ts object to a data frame
DecDate<- data.frame(Frcst = c(AP.predict), DecTime = c(time(AP.predict)))
head(DecDate)

#Convert decimal date to date format
DatesFormtd<- format(date_decimal(DecDate$DecTime), "%Y-%m-%d")
head(DatesFormtd)

#Bind calander date to data frame
Forecast<- cbind(DecDate, DatesFormtd)
head(Forecast)
plot(Forecast$DecTime, Forecast$Frcst, type="l", lty=1, lwd=2)

	#END TABLE OF FORECAST AND PLOT
####################################################
##10.3.2. ARIMA
####################################################
	#BEGIN ARIMA FORECAST

#Forecast with ARIMA model
d.arima<- auto.arima(AirPassengers)
d.forecast <- forecast(d.arima, level = c(95), h = 4*12)
autoplot(d.forecast)
str(d.arima)

	#END ARIMA FORECAST
####################################################	
	#BEGIN BUILD UPPER/LOWER/MEAN TABLE & PLOT

#Converts ts object to data frame
ArimaDf<- data.frame(Mean = c(d.forecast$mean), 
	Upper = c(d.forecast$upper), 
	Lower = c(d.forecast$lower), 
	DecTime = c(time(d.forecast$mean)))
head(ArimaDf)

#Convert decimal date to date format
DatesFormtdArima<- format(date_decimal(ArimaDf$DecTime), "%Y-%m-%d")
head(DatesFormtdArima)

ForecastArima<- cbind(ArimaDf, DateObj=DatesFormtdArima)
head(ForecastArima)

#Plot the forecast with upper and lower confidence interval
layout(1:1)
plot(ForecastArima$DecTime, ForecastArima$Mean, type="l", lty=1, lwd=2, ylim=c(400,900))
lines(ForecastArima$DecTime, ForecastArima$Upper, type="l", lty=2, lwd=2, col="red")
lines(ForecastArima$DecTime, ForecastArima$Lower, type="l", lty=2, lwd=2, col="red")

	#END BUILD UPPER/LOWER/MEAN TABLE & PLOT
####################################################
##10.3.3. TEST AS DETERMINISTIC OR STOCHASTIC
####################################################
	#BEGIN RANDOM

#Import dataset
StochSeries<- read.xlsx(navPath,
	sheetName="tblStochSeries", header=TRUE) 
head(StochSeries)

#Random series
layout(1:2)
RandomTs<- ts(StochSeries$Random, start = c(2007,1), freq=12)
plot(RandomTs, type ="l", lwd=2)
acf(RandomTs, lwd=2)

	#END RANDOM
####################################################
	#BEGIN RANDOM WALK

#Random walk series
layout(1:3)
RandomWalkTs<- ts(StochSeries$RandomWalk, start = c(2007,1), freq=12)
plot(RandomWalkTs, type ="l", lwd=2)
acf(RandomWalkTs, lwd=2)
acf(diff(RandomWalkTs), lwd=2)

#Compare to known deterministic series
layout(1:3)
plot(AP, lwd=2)
acf(AP, lwd=2)
acf(diff(AP, lwd=2))

	#END RANDOM WALK
####################################################
	#BEGIN DRIFT IN WALK

#Check for drift in Random walk
RandomWalkTsDrft<- diff(RandomWalkTs)
mean(diff(RandomWalkTsDrft))
mean(diff(RandomWalkTsDrft))+ 
	c(-1.96, 1.96)*sd(diff(RandomWalkTsDrft))/
	sqrt(length(RandomWalkTsDrft))

	#END DRIFT IN WALK
	#END SCRIPT