#install and load the forecast package
install.packages("forecast")
library(forecast)

#path is where the file is located , e.g. “/Users/tylermunger/Documents/tim245/”
time_series_data <- read.csv("$PATH$/time_series_data.csv")

time_series <- ts(time_series_data$value, start=c(2002, 1), end=c(2006, 4), frequency=4)

#plot the data, assuming time_series is a ts object
plot(time_series ,col="black", lwd=3)
plot(decompose(time_series))

#create the ARIMA models (see lecture notes)
arima_fit <- Arima(time_series, order= c(1,0,1), seasonal = c(1,1,0))
auto_arima_fit <- auto.arima(time_series, stepwise=FALSE, seasonal=TRUE, approximation=FALSE, D=1)

#create the exponential smoothing model (either holt or winters)
ets_fit <- ets(time_series)

#We can force a seasonal ets, i.e. winters, by explictly specifying the model type as shown below
ets_fit <- ets(time_series, model="ZZM")

#plot the models. red is the hard-coded arima model, blue is the auto-fit arima model, and green is the exponential smoothing model
plot(arima_fit$x,col="black", lwd=3)
lines(fitted(arima_fit),col="red")
lines(fitted(auto_arima_fit),col="blue")
lines(fitted(ets_fit),col="green")

#create the acf plots to determine the level of differencing
acf(time_series)

#create the differenced time series to account for the seasonal pattern
diff_time_series <- diff(time_series, 4)

#create the acf and pacf plots to determine the order of the AR and MA terms
acf(diff_time_series)
pacf(diff_time_series)

#error analysis of the models
summary(arima_fit)