Hwk 1 - Chapters 5 & 6 - Predictive Analytics

#5.10 Exercises

#assign daily 20 set
daily20 <- head(elecdaily,20)
frequency(daily20)

##.  electricity demand for Victoria, Australia, during 2014 is contained in elecdaily. The data for the first 20 days can be obtained as follows.
#1.a.Plot the data and find the regression model for Demand with temperature as an explanatory variable. Why is there a positive relationship? 
autoplot(daily20, facets=TRUE) + ggtitle("Time Series of Demand / Work Day / Temp")

daily20 %>%
  as.data.frame() %>%
  ggplot(aes(x=Temperature, y=Demand)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE)

#model
daily_elec_mod <- tslm(Demand ~ Temperature, data=daily20)
summary(daily_elec_mod)

#There is a positive relationship between temperature and demand because as hotter temperature lead to increased use of air cooling systems. It's also possible that at extreme levels, more people stay indoors.
#In fact as the half-hourly temp increases by 1 degrees, the demand for electricity likewise increases by 6.7gw

#1.b Produce a residual plot. Is the model adequate? Are there any outliers or influential observations?
checkresiduals(daily_elec_mod)

res_elec <- residuals(daily_elec_mod)
summary(res_elec)

#The good news is that the residuals appear to be white noise and not correlate as none cross the ACF threshold 
#It's also true that the mean of the residuals is essentially 0, so the two big criteria related to the mean of residuals and correlation are good
#The residuals however due show inconsistent variation around the mean with a left skewed, bimodal distribution
#The P-value from the Breush-Godfrey test also indicates that there's no reason to believe the residual are correlated any differently from white noise

#1.c Use the model to forecast the electricity demand that you would expect for the next day if the maximum temperature was
#15 degrees and compare with the forecast if the max was 35 degrees. Do you believe these forecasts?
forecast(daily_elec_mod, newdata = data.frame(Temperature = c(15, 35)))

#Both point forecasts [T(15) = 140.57 & T(35) = 275] trend as would be expected, meaning lower gw usage at lower temps
#Both figures would likely plot above the fitted line, which makes me wonder if some outlier influence is pulling our regression model lower

#1.d Give prediction intervals for your forecasts.
#For the 15 degree forecast, would expect the mean result of 95% of samples to fall within the interval of 108.6 and 190.9 degrees
#For the 35 degree forecast, would expect the mean result of 95% of samples to fall within the interval of 275.7 and 323.9 degrees

#1.e Plot Demand vs Temperature for all of the available data in elecdaily. What does this say about your model?
elecdaily %>%
  as.data.frame() %>%
  ggplot(aes(x=Temperature, y=Demand)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) + ggtitle("All Data Temp vs Demand")

#The plot of the entire data set as compared to the first 20 point subset indicates to me that our forecast has too high of a positive slope
#I think my model underestimates energy consumption as temperature drops when demand increases for heating 
#The forecast for 35 degrees still seems reasonable, while the 15 degree forecast is way understated.

#2. Data set mens400 contains the winning times (in seconds) for the men’s 400 meters final in each Olympic Games from 1896 to 2016.
#2a. Plot the winning time against the year. Describe the main features of the plot.
autoplot(mens400) + xlab("Year") + ylab("Time (s)") + ggtitle("Olympic Winning Times - Men's 400m")

#The main features of the plot are 1) a clear negative (or downward) trend from 1900 to 2016 in terms of the winning time to complete the race,
#2) Noticeable gaps in the results during WW1 & WW2, and 3)No seasonality

#2b. Fit a regression line to the data. Obviously the winning times have been decreasing, but at what average rate per year?
m400 = ts(mens400)
time400 = time(mens400)
m400_mod = tslm(formula = m400~time400, data = mens400)
autoplot(mens400, series = "Data") + autolayer(fitted(m400_mod), series = "Fitted") + 
  ggtitle("M400 Olympic Winning Times vs Regression Fit") + xlab("Year") + ylab("Time (s)")

summary(m400_mod)

#The intercept of 172.5 is not very meaningful in this instance as it assumes the trend is linear all the way back to year 0
#Essentially the intercept says that at the start of first millenium, the predicted 400m time was 172, again not meaningful
#However, the intercept tells us that each year sees a decrease in the winning 400m time by 0.06seconds, or
#Each 4 year olympic cycle, we would predict a ~0.25 second decline in the winning 400m time

#2c. Plot the residuals against the year. What does this indicate about the suitability of the fitted line?
checkresiduals(m400_mod)

summary(residuals(m400_mod))

#The plotted residuals indicates that though our R squared value is high, there remains significant predictor factors unexplained
#I have concerns that our fitted line has too much variance, despite passing the main criteria of no autocorrelation and mean 0 for residuals
#Autocorrelation within the residuals does not seem to be an issue as the Breusch-Godfrey test and ACF plot show as white noise
#The residuals show significant variance, notable at the start of the 20th century with a non-normal distribution and noteceable outlier

#2d. Predict the winning time for the men’s 400 meters final in the 2020 Olympics. Give a prediction interval for your forecasts. What assumptions have you made in these calculations?
forecast(m400_mod, newdata = data.frame(time400 = c(2020)))

#The predicted 400m winning time at the 2020 olympics is 42.04 seconds with 95% confidence interval between 39.5 and 44.5 seconds
#The actual winning time in Tokyo was 43/85 seconds, which would have been outside the 80% confidence interval

#3. Type easter(ausbeer) and interpret what you see.
easter(ausbeer)

#The data set contains quarterly columns (Q1:Q4) for row years from (1956:2010)
#The values in the matrix fall between 0 and 1
#These appear to represent dummy variables that indicate what quarter in which Easter is occurring as an explanatory factor for Austrailian beer sales
#I suspect that the fractional values represent when Holy Week is split between quarters

#4. Elasticity coefficient
#I'm not sure I fully understand the ask on this question
#First making y a function of x and showing what B1 equals
#Log(y) = B0 + B1Log(x) +e <--Starting equation
#Log(y) - B0 - e = B1Log(x)  <--Start to rearrange around B1Log(x)
#(Log(y)-B0 - e)/Log(x) = B1  <- Isolate elasticity coefficient
#B1 =  (log(y) - B0 - e)/log(x)
#Any percent change in preditor variable x would be equaled by the perecent change in forecast variable y

#5. The data set fancy concerns the monthly sales figures of a shop which opened in January 1987 and sells gifts, souvenirs, and novelties. The shop is situated on the wharf at a beach resort town in Queensland, Australia. The sales volume varies with the seasonal population of tourists. There is a large influx of visitors to the town at Christmas and for the local surfing festival, held every March since 1988. Over time, the shop has expanded its premises, range of products, and staff.
#5a.
head(fancy)

autoplot(fancy) + ylab("Sales volume") + xlab("Month, Years") + ggtitle("Monthly Sales Figure Fancy Shop")

ggsubseriesplot(fancy, main = NULL)

#There is clear seasonality at play in the data with annual spikes - towards the latter part of year (Q4 months)
#The seasonal spikes are also trending up (or increasing) as time advances
#March has a surprising peak relative to the rest of Q1 & Q2, which may be explained some extraordinary local event

#5b. Log transformations are useful when the data shows increases or decreases with the level of the series because these sets are difficult to fit a model
#The goal of the log transformation is to level out the seasonal variation over time
transformed_fancy <- BoxCox(fancy,0)
Original_fancy <- autoplot(fancy) + ylab("Sales volume") + xlab("Month, Years") + ggtitle("Original Fancy Shop")
Trans_fancy <- autoplot(transformed_fancy) +  xlab("Month, Years") + ylab("Log Sales Volume") + ggtitle("Transformed Fancy Shop")
grid.arrange(Original_fancy, Trans_fancy, ncol = 2)

#Learned about grid arrange scoping out how to facet multiple plots

#5c. Fit a regression model to the logarithms of these sales data with a linear trend and seasonal variables.
#I will add back in the trend and seasonality composite effects in the regression model
fancy_trans_mod <- tslm(transformed_fancy~trend+season)
summary(fancy_trans_mod)

#Did not add a specific dummy variable for surfing competition, creating forecasting issues in 5h. that I could not resolve

autoplot(transformed_fancy, series = "Transformed Data") + autolayer(fitted(fancy_trans_mod), series = "model") + ylab("Log Sales volume") + xlab("Month, Years") + ggtitle("Original Fancy Shop")

#5d. Plot the residuals against time and against the fitted values. Do these plots reveal any problems with the model?
autoplot(residuals(fancy_trans_mod)) + ylab("residuals")

fitted_res = cbind.data.frame(Fitted_values = fancy_trans_mod$fitted.values, Residuals = fancy_trans_mod$residuals)
ggplot(data = fitted_res, aes(x = Fitted_values, y = Residuals)) + geom_point()

#The analysis of the residuals seems to show clear indications of autocorrelation - it seems that seasonality has not been entirely neutralized,
#Heteroscedasticity is also in play as the variation between the residuals and fitted is not consistent

#5e.Do boxplots of the residuals for each month. Does this reveal any problems with the model?
#Month repeat for 7 years
residuals_month <- data.frame(fancy_trans_mod$residuals, month=rep(1:12,7))
boxplot(fancy_trans_mod$residuals ~ month, data=residuals_month)

#A better fitted model would show consistent variability in the residuals box plot, whereas March and the Fall show clear seasonal and cyclical effects

#5f. What do the values of the coefficients tell you about each variable?
summary(fancy_trans_mod)  

#The intercept indicates that the starting point in january is log sales volume of 7.6
#All subsequent coefficients are positive, which indicates that January is the lightest sales volume season
#March and Sept-Dec show strong sales.
#The model is trending in an increasing direction with each period

#5g.What does the Breusch-Godfrey test tell you about your model?
checkresiduals(fancy_trans_mod)

#The Breush-Godfrey test indicates that we can reject the hypothesis that our residuals are not correlated, or are just white noise
#This means that the model is not fully capturing the explainable effects

#5h. Regardless of your answers to the above questions, use your regression model to predict the monthly sales for 1994, 1995, and 1996. Produce prediction intervals for each of your forecasts.
fc_fancy <- forecast(fancy_trans_mod, h = 12*3)
autoplot((fc_fancy))+  xlab("Month, Years") + ylab("Log Sales Volume") + ggtitle("Fancy Shop Transformed Forecast 1994-1996")

#Shows prediction intervals
fc_fancy

#5i. Transform your predictions and intervals to obtain predictions and intervals for the raw data.
fc_raw_fancy <- exp(data.frame(fc_fancy))
ts_fc_raw_fancy <- ts(fc_raw_fancy, frequency = 12, start = 1994)
autoplot(ts_fc_raw_fancy) + xlab("Month, Years") + ylab("Sales Volume") + ggtitle("Fancy Shop Transformed Forecast 1994-1996")

#5j. How could you improve these predictions by modifying the model
#1) The first change I would make to my model is addressing the surf competition in March by adding dummy variables
#2) It seems like the model is faltering in terms of transformations appropriately reducing seasonal effects, so I would use a different log transformation

#6. The gasoline series consists of weekly data for supplies of US finished motor gasoline product, from 2 February 1991 to 20 January 2017. The units are in “million barrels per day.” Consider only the data to the end of 2004.
#End at 2005 to limit data to through 2004

#6a.
gas_2004 <- window(gasoline, end = 2005)
head(gas_2004)

#Weekly data, starts in 1991
str(gas_2004)

#Time series
autoplot(gas_2004) + xlab("Time (weekly") + ylab("Gas Volume (Million Barrels per Day)") + ggtitle("US Gas Supply on Weekly basis")

#Because the data shows weekly seasonality, Fourier modeling can be used to approximate periodic function
#Max fourier periods is 26 -> m/2
#For Fourier transformations, will try times of K= 4, 12, & 13

f.4_gas_2004 <- tslm(gas_2004 ~ trend + fourier(gas_2004,K = 4))
f.12_gas_2004 <- tslm(gas_2004 ~ trend + fourier(gas_2004,K = 12))
f.13_gas_2004 <- tslm(gas_2004 ~ trend + fourier(gas_2004,K = 13))

autoplot(gas_2004, series = "Data") + autolayer(fitted(f.4_gas_2004)) + autolayer(fitted(f.12_gas_2004)) + autolayer(fitted(f.13_gas_2004)) + xlab("Time (weekly)") + ylab("Gas Volume (Million Barrels per Day)") + ggtitle("US Gas Supply on Weekly basis")

#The Fourier harmonic models are much tighter in variance from the trend than the actual data

#6b. Select the appropriate number of Fourier terms to include by minimising the AICc or CV value.
#CV function displays 5 measures for prediction accuracy, with the goal of minimizing each, except adj R^2 (want higher)
CV(f.4_gas_2004)

CV(f.12_gas_2004)

CV(f.13_gas_2004)

#The harmonic series with a K value of 12, minimizes the cross-validation values and maximized adjR2
#I think f.12_gas_2004 is the best overall model, but if we wanted to minimize AIC, we'd select a lower K value model

#6c. Check the residuals of the final model using the checkresiduals() function. Even though the residuals fail the correlation tests, the results are probably not severe enough to make much difference to the forecasts and prediction intervals. (Note that the correlations are relatively small, even though they are significant.
checkresiduals(f.12_gas_2004)

#The residuals fit nicely around the mean 0 with relative normal distribution
#The residuals however do show some autocorrelation and fail the Breush-Gidfrey test indicating they are not simply whitenoise.

#6d. To forecast using harmonic regression, you will need to generate the future values of the Fourier terms. This can be done as follows.
#Provided -> fc <- forecast(fit, newdata=data.frame(fourier(x,K,h)))
fc_gas_2005 <- forecast(f.12_gas_2004, newdata = data.frame(fourier(gas_2004, K=12, h = 52)))
fc_gas_2005

#6e. Plot the forecasts along with the actual data for 2005. What do you find?
gas_2005_actuals <- window(gasoline, start = 2005, end = 2006)
autoplot(gas_2005_actuals, series = "Actual data") + autolayer(fc_gas_2005$mean, series = "Prediction Mean") + xlab("Time (weekly)") + ylab("Gas Volume (Million Barrels per Day)") + ggtitle("US Gas Supply 2005: Actual vs Prediction")                        

autoplot(gas_2005_actuals, series = "Actual data") + autolayer(fc_gas_2005, series = "Prediction Range")                        

autoplot(fc_gas_2005, series = "Prediction Range") + autolayer(gas_2005_actuals, series = "Actual data") + autolayer(fc_gas_2005$mean, series = "Prediction Mean")

#The prediction follows the 2005 actual trend but in a less extreme manner (highs and lows are smaller)
#The drop in actual gas volume in Q3 of 2005 is not captured within the 95% confidence interval for prediction

#7. Data set huron gives the water level of Lake Huron in feet from 1875 to 1972.

#7a. Plot the data and comment on its features.
#Huron lake level in feet by year
frequency(huron)

autoplot(huron) + xlab("Year") + ylab("Huron Depth (ft)") + ggtitle("Depth of Lake Huron by Year")

#The depth of the lake peaked just before 1880 and hit a nadir in ~1965
#The trend does seem to be decreasing with periodic bouncebacks
#I would expect there to be cyclical movement related to periods of draught or heavy precipitation

huron_45 <- window(huron, start = 1945, end = 1950)
autoplot(huron_45)

#Minimal seasonality when looking at a random 5 year period

#7b. Fit a linear regression and compare this to a piecewise linear trend model with a knot at 1915.
#leaving out seasonality as it does not seem to play a role

#linaer model
lin_mod_huron <- tslm(huron~trend)

#Piecewise model
Huron_T <- time(huron,na.action = na.omit) 
TS_PW <- ts(pmax(0,Huron_T-1915), start=1875)
Huron_pw_mod <- tslm(huron~Huron_T+TS_PW)

#Model Summaries
summary(lin_mod_huron)

summary(Huron_pw_mod)

checkresiduals(lin_mod_huron)

checkresiduals(Huron_pw_mod)

autoplot(huron, series = "Data") + autolayer(fitted(lin_mod_huron), series = "Linear Reg") + autolayer(fitted(Huron_pw_mod), series = "Piece-wise Mod") + xlab("Year") + ylab("Huron Depth (ft)") + ggtitle("Depth of Lake Huron by Year")

#The adj R^2 value is improved by the piecewise function as it follows the downward trend to 1915 where it flattens the trend

#7c. Generate forecasts from these two models for the period up to 1980 and comment on these.
fc_lm <- forecast(lin_mod_huron, h = 8)
fc_lm

autoplot(huron, series = "Data") + autolayer(fc_lm, series = "Lin Model") + autolayer(fitted(fc_lm),series = "Lin Model") + xlab("Year") + ylab("Huron Depth (ft)") + ggtitle("Linear Forecast 1972-1980")

#The piecewise forecast requires some manipulation of the new data period

fc_time <- Huron_T[length(Huron_T)] + seq(8)
fc_ts <- TS_PW[length(TS_PW)] + seq(8)
fc_pw <- forecast(Huron_pw_mod, newdata = as.data.frame(cbind(Huron_T = fc_time, TS_PW = fc_ts)))
autoplot(huron, series = "Data") + autolayer(fc_pw, series = "PW Model") + autolayer(fitted(fc_pw),series = "PW Model") + xlab("Year") + ylab("Huron Depth (ft)") + ggtitle("PW Forecast 1972-1980")

#The linear model more effectively captures the most recent actual data in its confidence interval
#But the PW model more accurately reflect the leveled off trending of the data

#8 Sorry- I did not have capacity to delve into the matrices at this time

#CHAPTER 6 Exercises

#1 Show that a 3 × 5 MA is equivalent to a 7-term weighted moving average with weights of 0.067, 0.133, 0.200, 0.200, 0.200, 0.133, and 0.067.
#The 3 x 5 MA is equivalent to 1/3[1/5(y1 + Y2 + y3 + y4 + y5) + 1/5(Y2 + y3 + y4 + y5 + y6) + 1/5(y3 + y4 + y5 + y6+y7)]
#When factored out -> (1/15)y1 + (2/15)y2 + (3/15 = 1/3)y3 + (3/15 = 1/3)y4 + (3/15 = 1/3)y5 + (2/15)y6 + 1/15y7
y1 = 1/15
y2 = 2/15
y3 = 3/15
y4 = 3/15
y5 = 3/15
y6 = 2/15
y7 = 1/15
MA_6.1 <- c(y1, y2, y3, y4, y5, y6, y7)
MA_6.1

#The resulting 3 x 5 MA match the 7-term weighted moving average

#2 The plastics data set consists of the monthly sales (in thousands) of product A for a plastics manufacturer for five years.

#2a. Plot the time series of sales of product A. Can you identify seasonal fluctuations and/or a trend-cycle?
#?plastics
#Monthly sales of product A for a plastics manufacturer
frequency(plastics)

#Frequency is months
autoplot(plastics)

#An overall upward trend that shows increasing sales volume from year 1 - 5
ggseasonplot(plastics)

plastics_month <- data.frame(plastics, month=rep(1:12,5))
boxplot(plastics ~ month, data=plastics_month)

#A seasonal trend is also clearly present across all 5 years, an initial decline from January to Feb followed by increasing sales that peak during the Summer and begins to decline in Aug - Dec

#2b. Use a classical multiplicative decomposition to calculate the trend-cycle and seasonal indices.
plastics_multdecomp <- decompose(plastics, type = "multiplicative")
autoplot(plastics_multdecomp, Facet = TRUE)

#2c.Do the results support the graphical interpretation from part a?
#The plot confirms the analysis from 2.a. 
#The trend is a clear increase from year 1->5 with some later leveling
#Seasonal trends are consistent across years with clear peaks and troughs

#2d.Compute and plot the seasonally adjusted data. 
autoplot(seasadj(plastics_multdecomp)) + ylab("Plastic Sales Volume") + xlab("Months, Years")

#The removal of the seasonal effect leaves just the upwards trend in sales

#2e. Change one observation to be an outlier (e.g., add 500 to one observation), and recompute the seasonally adjusted data. What is the effect of the outlier?
plastics_outlier <- plastics
plastics_outlier[36] = plastics_outlier[36] + 750
autoplot(plastics_outlier, series = "Outlier") + autolayer(plastics, series = "Original")

plastics_multdecomp_outlier <- decompose(plastics_outlier, type = "multiplicative")
autoplot(plastics_multdecomp_outlier, Facet = TRUE) 

autoplot(seasadj(plastics_multdecomp_outlier)) + ylab("Plastic Sales Volume") + xlab("Months, Years")

#The outlier continues to impact the data and is not resolved by the seasonal adjustment

#2f. Does it make any difference if the outlier is near the end rather than in the middle of the time series?
#Adjustments for seasonal changes will not resolve the outlier regardless if it is the middle or near the end of the data set.

#3. Recall your retail time series data (from Exercise 3 in Section 2.10). Decompose the series using X11. Does it reveal any outliers, or unusual features that you had not noticed previously?
getwd()

retaildata <- read_xlsx("retail.xlsx", skip=1)
head(retaildata)

#My name starts with a T, so I am choosing the data column with a T in it
retaildata_ts <- ts(retaildata[,"A3349335T"], frequency=12, start=c(1982,4))
#Running some of the initial analysis I would have done from section 2
autoplot(retaildata_ts)

ggseasonplot(retaildata_ts)

#Clear increasing trend in retail sales by year
#Seasonal spike in more recent years in March - could this be related to a special event
#Seasonal increases during holiday season - November & December
retaildata_x11 <- seas(retaildata_ts, x11(""))
autoplot(retaildata_x11, facet = TRUE)

#The two key noticings from the x11 analysis
#1)The seasonal effect was larger before 1990
#2)Some large irregularities (maybe outliers) occurred around 1986 and 2009 that could impact modeling

#4. Figures 6.16 and 6.17 show the result of decomposing the number of persons in the civilian labour force in Australia each month from February 1978 to August 1995.

#4a.Write about 3–5 sentences describing the results of the decomposition. Pay particular attention to the scales of the graphs in making your interpretation.
#The Austrailian labor force data shows a clear positive trend so that the level of civilian employment is increasing over time. There was a noticeable outlier dip in the early 1990s which could be the result of a more general business disruption, like a market correction. The trend line confirms this analysis with more steady growth from 1980-1985 and more rapid expansion 1985-1990.  Seasonal employment increases seemed to happen during the June and July than dip in August, which may be the result of school break employment.

#4b. Is the recession of 1991/1992 visible in the estimated components?
#The remainder components during the years 1991-1992 indicate that actual employment levels were far below composite modeling, which indicates a recession did occur.

#5. This exercise uses the cangas data (monthly Canadian gas production in billions of cubic metres, January 1960 – February 2005).
#5a. Plot the data using autoplot(), ggsubseriesplot() and ggseasonplot() to look at the effect of the changing seasonality over time. What do you think is causing it to change so much?
autoplot(cangas)

ggseasonplot(cangas)

ggsubseriesplot(cangas, TRUE, TRUE)

#There's a general upwards trend that levels between 1975-1990 with larger increases from 1960-1975 & 1990-2005
#Regarding seasonality, colder months Sept-Feb show higher gas usage, which is likely the result of increased use for heating homes
#The general upward trend reflection increasing population and increasing national wealth to afford heating

#5b. Do an STL decomposition of the data. You will need to choose s.window to allow for the changing shape of the seasonal component.
cangas_window <- stl(cangas, s.window = 12)
autoplot(cangas_window)

#The trend in the stl decomp confirms the intial analysis that the trend flattens in the early 80s, while also revealing that the seasonal impact increased in the '80s

#5c. Compare the results with those obtained using SEATS and X11. How are they different?
cangas_x11<-seas(cangas, x11 = "")
autoplot(cangas_x11)

#The trend factor seems to be less smooth then with STL decomp, while the scale of seasonal impact narrower but more significant in different periods.  For example the 1960s are larger effects, while the 1980s are lessened.

#6 We will use the bricksq data (Australian quarterly clay brick production. 1956–1994) for this exercise.
#?bricksq -> Australian quarterly clay brick production: 1956–1994.

#6a. Use an STL decomposition to calculate the trend-cycle and seasonal indices. (Experiment with having fixed or changing seasonality.)
bricksq_stl_period <- stl(bricksq, s.window = "periodic", TRUE)
bricks_period_plot <- autoplot(bricksq_stl_period) + ggtitle("Periodic STL")
bricksq_stl_fixed <- stl(bricksq, s.window = 12)
bricks_fixed_plot <- autoplot(bricksq_stl_fixed) +ggtitle("Fixed STL (12)")
grid.arrange(bricks_period_plot,bricks_fixed_plot, ncol = 2)

#The fixed plot creates a less consistent seasonal effect

#6b. Compute and plot the seasonally adjusted data.
brick_mod <- seasadj(stl(bricksq, t.window = 12, s.window = "periodic", robust = TRUE))
brick_mod_plot <- autoplot(brick_mod) + ggtitle("Brick Model - Seasonal Adj")
brick_plot <- autoplot(bricksq) + ggtitle("Brick Data")
grid.arrange(brick_plot, brick_mod_plot, ncol = 2)

# Seasonal adj is smoother

#6c. Use a naïve method to produce forecasts of the seasonally adjusted data.
brick_mod_naive <- naive(brick_mod)
autoplot(brick_mod_naive)

#6d. Use stlf() to reseasonalise the results, giving forecasts for the original data.
#Adding random walk to stl decomp
bricksq_stlf <- stlf(bricksq, method = "naive")
stlf_mod <- forecast(bricksq_stlf)
autoplot(stlf_mod)

#6e/f Compare residuals from STL and random walk
checkresiduals(brick_mod_naive)

checkresiduals(stlf_mod)

#Both sets of residuals have a left tail distribution with inconsistent variance
#Naive method show lower variance, but marginally higher autocorrelation than random walk

#6g Compare forecasts from stlf() with those from snaive(), using a test set comprising the last 2 years of data. Which is better?
train_brick <- window(bricksq, end = c(1992,3))
brick_snaive_mod <- snaive(train_brick)
brick_stlf_mod <- stlf(train_brick)
brick_snaive_fc <- forecast(brick_snaive_mod, h=8)
brick_stlf_fc <- forecast(brick_stlf_mod, h=8)
autoplot(bricksq, series = "data") + autolayer(fitted(brick_snaive_fc), series = "SNaive") + autolayer(brick_snaive_fc, series = "SNaive") + ggtitle("SNaive vs Actuals")

autoplot(bricksq, series = "data") + autolayer(fitted(brick_stlf_fc), series = "STLF") + autolayer(brick_stlf_fc, series = "STLF")+ ggtitle("STLF vs Actuals")

accuracy(brick_snaive_fc, bricksq)

accuracy(brick_stlf_fc, bricksq)

checkresiduals(brick_snaive_mod)

checkresiduals(brick_stlf_mod)

#The random walk (STLF) decomp has significant advantages in terms of lower autocorrelation and smaller error factors

#7. Use stlf() to produce forecasts of the writing series with either method="naive" or method="rwdrift", whichever is most appropriate. Use the lambda argument if you think a Box-Cox transformation is required.
head(writing)

#Industry sales for printing and writing paper in thousands of Francs - 1963-1972
autoplot(writing)

ggseasonplot(writing)

#Strong seasonality, converges to almost nil every August
#Because of the strong seasonal effects, I will utilize the box-cox transformations
writing_stlf_naive <- stlf(writing, s.window = 12, method = "naive", lambda = BoxCox.lambda(writing))
writing_stlf_drift <- stlf(writing, s.window = 12, method = "rwdrift", lambda = BoxCox.lambda(writing))
autoplot(writing_stlf_naive) + ylab("Log Industry Sales (100s Francs") + xlab("years, months")

 autoplot(writing_stlf_drift) + ylab("Log Industry Sales (100s Francs") + xlab("years, months")

checkresiduals(writing_stlf_naive)

checkresiduals(writing_stlf_drift)

#Both the random walk with and without drift show very similar forecasts
#The residuals for the drift form look slightly more normal and the upward trend is slightly acentuated so I prefer it in this case

#8 Use stlf() to produce forecasts of the fancy series with either method="naive" or method="rwdrift", whichever is most appropriate. Use the lambda argument if you think a Box-Cox transformation is required.
#For question 7 I compared naive and drift, here I am going stick with the drift method and compare with and without a Box.cox transformation
#See chapter 5 question 5 for previous Fancy data analysis  
fancy_stlf_drift <- stlf(fancy, s.window = 12, method = "rwdrift")
fancy_stlf_drift_boxcox <- stlf(fancy, s.window = 12, method = "rwdrift", lambda = BoxCox.lambda(writing))
fsd_plot <- autoplot(fancy_stlf_drift) + ylab("Sales volume")
fsd_bc_plot <- autoplot(fancy_stlf_drift_boxcox) + ylab("Log Sales Volume")
grid.arrange(fsd_plot, fsd_bc_plot, ncol=2)

checkresiduals(fancy_stlf_drift)

checkresiduals(fancy_stlf_drift_boxcox)

#The boxcox transformation helps normalize for various outliers