Forecasting Assignment Week 2

## R Packages needed
library(readr)
library(readxl)
library(ggplot2)
library(fpp2)
library(psych)
library(lmtest)
library(seasonal)


#Chapter 5
#Q1
daily20 <- head(elecdaily,20)

#a.
autoplot(daily20, main="Electricity Demand for Victoria, Australia during 2014")
lm(Demand~Temperature, data=elecdaily)

Call:
lm(formula = Demand ~ Temperature, data = elecdaily)

Coefficients:
(Intercept)  Temperature  
   212.3856       0.4182  

summary(lm(Demand~Temperature, data=elecdaily))

Call:
lm(formula = Demand ~ Temperature, data = elecdaily)

Residuals:
    Min      1Q  Median      3Q     Max 
-55.474 -15.719  -0.336  15.767 117.184 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 212.3856     5.0080  42.409   <2e-16 ***
Temperature   0.4182     0.2263   1.848   0.0654 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.55 on 363 degrees of freedom
Multiple R-squared:  0.009322,  Adjusted R-squared:  0.006592 
F-statistic: 3.416 on 1 and 363 DF,  p-value: 0.0654

#Demand and temp have a low correlation. The p value of temp is considered a explantory variable of demand and has the value of .0654. temp shows a small increase throughout different points of the plot which causes high demand. the biggest increase we see in demand is when temp increases around "time 3"
#b.
plot(lm(Demand~Temperature, data=elecdaily))

#It looks like a majority of values are outside ofth the cooks distance line which means the values have little value. looking at observation 15-17 you see that they are valuable to the regression results and shold be looked at more. for this plot it would be best to remove outliers

#c

DemTSLM <- tslm(Demand ~ Temperature, data=daily20)
tempRange <- data.frame(Temperature=c(15,35))
forecastDemTemp <- forecast(DemTSLM, tempRange)
forecastDemTemp

#the forecast seems reasonable especialy when compared to the deamnd values  from the training set. similat results can be found if you ran this forecast with the training data set.

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

fit <- tslm(Demand ~ Temperature, data=daily20)
checkresiduals(fit)

    Breusch-Godfrey test for serial correlation of order up to
    5

data:  Residuals from Linear regression model
LM test = 3.8079, df = 5, p-value = 0.5774

forecast(fit, newdata=data.frame(Temperature=c(15,35)))

#these prediction values are reasonable estimates when compared to the training set.

#e
plot(Demand ~ Temperature, data=elecdaily, main="Demand vs. Temperature of Elecdaily dataset")

#the plot above resembles the shape of the daily20 dataframe. The two sets of data share outliers from the residuals section. This proves that the subset of daily20 is a accurate representation of the full set of observations.

#q2
autoplot(mens400, main="Winning times for the Men's 400 meter Olympic Games Final", ylab="Time (Seconds)", xlab="Year")

#there is definately a inverse relationship between race time per second and year. the time to win the 400 m dash had decreased consistanly over time. there are several spots on the graph that are blank due to NA values being in the data.

#b.
#change dataset to time series format
time400 <- time(mens400)
tslm(mens400 ~ time400, data=mens400)

Call:
tslm(formula = mens400 ~ time400, data = mens400)

Coefficients:
(Intercept)      time400  
  172.48148     -0.06457  

#winning time has decreased by .06457 each year.
#add the regression line
autoplot(mens400, main="Winning times for the Men's 400 meter Olympic Games Final", ylab="Time (Seconds)", xlab="Year") + geom_smooth(method="lm", se=FALSE)

#c.
mens400RG <- lm(mens400 ~ time400, data=mens400)
plot(mens400RG)

#the residuals vs. fitted graph show that the slow is decreasing and the winning persons times decrease at a decreasing rate each year. this shows that the fitted line is accurate for shortterm. This doesnt how ever help forecast future winning times.

#d.

TimeMens400 <- time(mens400)
mens400LM <- lm(mens400 ~ TimeMens400, data=mens400, na.action=na.exclude)

#e. remove the 3 na values
Mens400F2022 <- forecast(mens400LM, newdata=data.frame(TimeMens400=2020))
Mens400F2022

#q3.
easter(ausbeer)

     Qtr1 Qtr2 Qtr3 Qtr4
1956 0.67 0.33 0.00 0.00
1957 0.00 1.00 0.00 0.00
1958 0.00 1.00 0.00 0.00
1959 1.00 0.00 0.00 0.00
1960 0.00 1.00 0.00 0.00
1961 0.33 0.67 0.00 0.00
1962 0.00 1.00 0.00 0.00
1963 0.00 1.00 0.00 0.00
1964 1.00 0.00 0.00 0.00
1965 0.00 1.00 0.00 0.00
1966 0.00 1.00 0.00 0.00
1967 1.00 0.00 0.00 0.00
1968 0.00 1.00 0.00 0.00
1969 0.00 1.00 0.00 0.00
1970 1.00 0.00 0.00 0.00
1971 0.00 1.00 0.00 0.00
1972 0.33 0.67 0.00 0.00
1973 0.00 1.00 0.00 0.00
1974 0.00 1.00 0.00 0.00
1975 1.00 0.00 0.00 0.00
1976 0.00 1.00 0.00 0.00
1977 0.00 1.00 0.00 0.00
1978 1.00 0.00 0.00 0.00
1979 0.00 1.00 0.00 0.00
1980 0.00 1.00 0.00 0.00
1981 0.00 1.00 0.00 0.00
1982 0.00 1.00 0.00 0.00
1983 0.00 1.00 0.00 0.00
1984 0.00 1.00 0.00 0.00
1985 0.00 1.00 0.00 0.00
1986 1.00 0.00 0.00 0.00
1987 0.00 1.00 0.00 0.00
1988 0.00 1.00 0.00 0.00
1989 1.00 0.00 0.00 0.00
1990 0.00 1.00 0.00 0.00
1991 1.00 0.00 0.00 0.00
1992 0.00 1.00 0.00 0.00
1993 0.00 1.00 0.00 0.00
1994 0.00 1.00 0.00 0.00
1995 0.00 1.00 0.00 0.00
1996 0.00 1.00 0.00 0.00
1997 1.00 0.00 0.00 0.00
1998 0.00 1.00 0.00 0.00
1999 0.00 1.00 0.00 0.00
2000 0.00 1.00 0.00 0.00
2001 0.00 1.00 0.00 0.00
2002 1.00 0.00 0.00 0.00
2003 0.00 1.00 0.00 0.00
2004 0.00 1.00 0.00 0.00
2005 1.00 0.00 0.00 0.00
2006 0.00 1.00 0.00 0.00
2007 0.00 1.00 0.00 0.00
2008 1.00 0.00 0.00 0.00
2009 0.00 1.00 0.00 0.00
2010 0.00 1.00          

# I see a data set that has the quarterly values from 1956 to 2010, with the ending in q2 of 2010. Each value for the quarters seem to be between 0 and 1. Could it signify a percentage of probability of a specific event happening?

#Q4
#log(y) = β0 + β1log(x) + e
#β1log(x) = log(y) - β0 - e
#β1 = (log(y) - β0 - e) / log(x)              

#q5
autoplot(fancy, main="Sales Volume over Time", ylab="Sales Volume")
seasonplot(fancy, main="Seasonal Sales Volume over Time", ylab="Sales Volume")

#the graph shows seasonal trends in the sales. In each year around Nov and Dec you can see the spike in sales due to the holiday season. we also can see that in 1992 and 1993 something happened that caused extremely high growth and expansion.

#b.
#Patterns in data sets that are plotted can be tricky and may show heteroscedasticity errors. Meaning that the variance of the residuals may not be constant and actually could lead to more issues as you use the data to run more advanced regression.The transformation of using a logarith can help combat this issue.

#c.
## Creating a "surfing festival" dummy variable
time <- time(fancy)

surfingFestival <- c()
for(i in 1:length(time)){
  month <- round(12*(time[i] - floor(time[i]))) + 1
  year <- floor(time[i])
  if(year >= 1988 & month == 3) ## started March, 1998
  {
    surfingFestival[i] <- 1 ## Binary - 1 for if surfing festival is in town
  } else {
    surfingFestival[i] <- 0 ## 0 if the surfing festival is not in town
  }
}

## Regression formula and output
## Using BoxCox power transformation prior to regression formula
fancyTrans <- (BoxCox(fancy, 0))
fancyReg <- tslm(fancyTrans ~ trend + season + surfingFestival)
summary(fancyReg)

Call:
tslm(formula = fancyTrans ~ trend + season + surfingFestival)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.33673 -0.12757  0.00257  0.10911  0.37671 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     7.6196670  0.0742471 102.626  < 2e-16 ***
trend           0.0220198  0.0008268  26.634  < 2e-16 ***
season2         0.2514168  0.0956790   2.628 0.010555 *  
season3         0.2660828  0.1934044   1.376 0.173275    
season4         0.3840535  0.0957075   4.013 0.000148 ***
season5         0.4094870  0.0957325   4.277 5.88e-05 ***
season6         0.4488283  0.0957647   4.687 1.33e-05 ***
season7         0.6104545  0.0958039   6.372 1.71e-08 ***
season8         0.5879644  0.0958503   6.134 4.53e-08 ***
season9         0.6693299  0.0959037   6.979 1.36e-09 ***
season10        0.7473919  0.0959643   7.788 4.48e-11 ***
season11        1.2067479  0.0960319  12.566  < 2e-16 ***
season12        1.9622412  0.0961066  20.417  < 2e-16 ***
surfingFestival 0.5015151  0.1964273   2.553 0.012856 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.179 on 70 degrees of freedom
Multiple R-squared:  0.9567,    Adjusted R-squared:  0.9487 
F-statistic:   119 on 13 and 70 DF,  p-value: < 2.2e-16

#d.
autoplot(fancyReg$residuals, main="Residuals plot of fancyReg Resgression", ylab="Residuals")
#looking at the plot the seasonality pattern exist still after the transformation. this is an issue for autocorrelation

#e
boxplot(fancyReg$residuals, main="fancyReg Regression Residuals")

#f.
fancyReg$coefficients

    (Intercept)           trend         season2         season3 
     7.61966701      0.02201983      0.25141682      0.26608280 
        season4         season5         season6         season7 
     0.38405351      0.40948697      0.44882828      0.61045453 
        season8         season9        season10        season11 
     0.58796443      0.66932985      0.74739195      1.20674790 
       season12 surfingFestival 
     1.96224123      0.50151509 

#the surfing festival coefficent has a value of .5015 which shows a relative increase in sales during march. The trend coefficent is fairly low at .0220 but does show a small increase in sales. we definately see se a large increase in the coefficent around season 10 which is probably due to the result of the shop expansion during the time.
library(lmtest)
bgtest(fancyReg)

    Breusch-Godfrey test for serial correlation of order up to
    1

data:  fancyReg
LM test = 25.031, df = 1, p-value = 5.642e-07

#the output from the BG test tells a story about serial correlation in the regression model. The p value is extremely low which tells us the null hypothesis can be rejected. That does however mean heteroscedasticity exists in my model.

#h. 
timeRangeFancy <- ts(data=time, start=1994, end=c(1996,12), frequency=12)
forecast(fancyReg, timeRangeFancy)

newdata column names not specified, defaulting to first variable required.

#i.
fancyTrans <- (BoxCox(fancy, 0))
forecast(fancyTrans)

#j.
#the lambda function could start to improve these functions. the function BoxCox.Lambda(fancy) would be a good way to transform a regression model. A better transformation can caputre more exponential growth better.

#q6.
#a. I could not find the dataset easily according to documentation from the text is actually is named euretail. which is what i will end up using.
gasoline2014 <- window(euretail, end=2005)
autoplot(gasoline2014, main="Quarterly Retail Trade Index in the Euro area", xlab="Year", ylab="Thousand barrels per day")

gasoline %>% tbats() %>% forecast() %>%
  autoplot() + xlab("Year") + ylab("Thousand barrels per day")

#b.
gasolineCheck <- tslm(gasoline2014 ~ trend)
gasolineCheck

Call:
tslm(formula = gasoline2014 ~ trend)

Coefficients:
(Intercept)        trend  
     88.610        0.302  

#c.
checkresiduals(gasolineCheck)

    Breusch-Godfrey test for serial correlation of order up to
    7

data:  Residuals from Linear regression model
LM test = 21.937, df = 7, p-value = 0.002604

#d.
gasolineF <- forecast(gasoline2014)
gasolineF

#e.
autoplot(gasolineF)

#the range of the forecast is pretty small which will lead to a more accurate forecast. using the 95% internal makes a much smaller / tighter prediction range.

#q7
autoplot(huron, main="Mean July Average Water Surface Elevation", ylab="Feet", xlab="Year")

# the plot shows a negative trend in water level between 1880 and 1986. This data also has seasonaility. The lowest level was around 1960 and the highest being around 1880.

#b

tslm(huron ~ trend)

Call:
tslm(formula = huron ~ trend)

Coefficients:
(Intercept)        trend  
    10.2020      -0.0242  

huronRG1 <- tslm(huron ~ trend)

tHuron <- time(huron)
tHuron15 <- 1915

tHuronSection <- ts(pmax(0,tHuron-tHuron15), start=1875)

#c

## Forecast 1
huronRG1 <- tslm(huron ~ trend)
forecast(huronRG1, h=8)

## Forecast 2
huronRG2 <- tslm(huron ~ tHuronSection)
forecast(tHuronSection, h=8)

# the two forecast swho some differences but both have different formatting aswell. Overall i see a increasing trend in both.



#chapter 6
#q1

## Centered moving averages can be smoothed by another moving average. This creates a double moving average. In the case of a 3x5 moving average, this signifies a 3 moving average of a 5 moving average. 

## Weights = c(0.067, 0.133, 0.200, 0.200, 0.200, 0.133, 0.067)

## 3x5 MA = [((Y1 + Y2 + Y3 + Y4 + Y5)/5) + ((Y2 + Y3 + Y4 + Y5 + Y6)/5) + ((Y3 + Y4 + Y5 + Y6 + Y7)/5)] / 3
## Plugging in these values proves that the 3x5 moving average is equal to a 7-term weighted moving average


#q2
#a.
autoplot(plastics, main="Monthly sales of Product A", xlab="Year", ylab="Sales (thousands")

#the plot shows seasonality of sale in the beginning and end parts of the year. sales then to be lower in the beginning of the year and higher towards the end.

#b
plastics %>% decompose(type="multiplicative") %>%
  autoplot() + xlab("Year") +
  ggtitle("Monthly sales of Product A")

#c.

#the results definately do support the statements from part a. it shows both the seasonality and the slight upward positive trend.

#d
autoplot(plastics, main="Monthly sales of Product A", ylab="Sales (Thousands)", xlab="Year") + autolayer(snaive(plastics, h=30), series="Seasonal Naïve", PI=FALSE) + autolayer(naive(plastics, h=30), series="Naïve", PI=FALSE) + autolayer(rwf(plastics, h=30), series="Drift", PI=FALSE)

#e
plasticsAdd <- plastics
plasticsAdd[15] <- plasticsAdd[15] + 500

autoplot(plasticsAdd, main="Monthly sales of Product A (Modified Value 15)", ylab="Sales (Thousands)", xlab="Year") + autolayer(snaive(plasticsAdd, h=30), series="Seasonal Naïve", PI=FALSE) + autolayer(naive(plasticsAdd, h=30), series="Naïve", PI=FALSE) + autolayer(rwf(plasticsAdd, h=30), series="Drift", PI=FALSE)

#the values forecasted stay stable with only the addition of 500 units to value 15. this supports the accurate nature  of the seasonally adjusted data prediction.


#f.
## Add 500 to end Value
plasticsAddEnd <- plastics
plasticsAddEnd[40] <- plasticsAddEnd[40] + 500

autoplot(plasticsAddEnd, main="Monthly sales of Product A (Modified Value 40)", ylab="Sales (Thousands)", xlab="Year") + autolayer(snaive(plasticsAddEnd, h=30), series="Seasonal Naïve", PI=FALSE) + autolayer(naive(plasticsAddEnd, h=30), series="Naïve", PI=FALSE) + autolayer(rwf(plasticsAddEnd, h=30), series="Drift", PI=FALSE)
# adding an outlier probably shifted the end value up slightly but not in a significant manner. overall outliers should always be evaluated for validity when using multiplicative decomposition.

#3.
library(readxl)
retail <- read_excel("R/retail.xlsx",skip=1)
## Creating a time series of the imported data
retailTS <- ts(retail[,"A3349873A"], frequency=12, start=c(1982, 3))
#plot
autoplot(retailTS, main="Monthly Austrailian Retail Sales", ylab="Sales", xlab="Year")

## X11 Decomposition
library(seasonal)
retailx11 <- seas(retailTS, x11="")

autoplot(retailx11, main="X11 Decomposition on Monthly Austrailian Retail Sales", xlab="Year")

#outliers are shown in the bottom remainder plot. Major outliers are 1989,'94,2001,2012.

#q4
#a

#Decomposition shows the results of the # of people in the civilian labour force from Austrialia each month from Feb '78 to Aug '95. Right off the bat we see a strong postive trend of the number of workers during this time period. The seasonality chart does show the growth is cyclical and has strong seasonality. There also are some major outliers that are in the data around '92. I think more research would be needed to explain the increase in labourers during this period. the second chart shows a very interesting look to the seasonal part of the months. We see a large increase in july and huge decrease in March.

#b.

#Above i talked about the large outliers, they would actually perfectly explain the results we see. The remainder plot does a great job at showing this.


#q5
#a
autoplot(cangas, main="Monthly Canadian Gas Production", ylab="Gas Production (Billions)", xlab="Year")

ggsubseriesplot(cangas, main="Monthly Canadian Gas Production", ylab="Gas Production (Billions)")
ggseasonplot(cangas, main="Seasonal Plot: Monthly Canadian Gas Production", ylab="Gas Production (Billions)")

#The seasonality probably is due to the fuel prices and demand during specific seasons. For example fuel consumption goes up in the summer since people travel more, which would lead to higher oil production. 


#b.
cangas %>%
  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot()


#c.
# The results are pretty similar but some there are some differences. the seasonality graph shows consistent sesaonality throught the years. The remainder plot provides some insight into the outlires that can be seen between '80 and '90. smaller outliers also appear between 1995 and 2005.


#q6
#a
bricksq %>%
  stl(t.window=26, s.window="periodic", robust=TRUE) %>%
  autoplot()

#b
autoplot(bricksq, main="Monthly sales of Product A", ylab="Sales (Thousands)", xlab="Year") + autolayer(snaive(bricksq, h=30), series="Seasonal Naïve", PI=FALSE)

#c
autoplot(bricksq, main="Australian quarterly clay brick production. 1956–1994", ylab="Brick Production", xlab="Year")+ autolayer(naive(bricksq, h=30), series="Naïve", PI=FALSE)

#d
brickF1 <- stlf(bricksq)
brickF1

autoplot(brickF1)


#e
checkresiduals((brickF1))

    Ljung-Box test

data:  Residuals from STL +  ETS(M,N,N)
Q* = 41.128, df = 6, p-value = 2.733e-07

Model df: 2.   Total lags used: 8

#f
brickF1R <- stlf(bricksq, robust=TRUE)
brickF1R

autoplot(brickF1R)
checkresiduals(brickF1R)

    Ljung-Box test

data:  Residuals from STL +  ETS(M,N,N)
Q* = 28.163, df = 6, p-value = 8.755e-05

Model df: 2.   Total lags used: 8

#those functions seem to have reduced normality but the residuals still have autocorrelation issues.


#g
trainBrick <- subset(bricksq, end=length(bricksq) - 9)
testBrick <- subset(bricksq, start=length(bricksq) - 8)

sBrick <- snaive(trainBrick)
stlfBrick <- stlf(trainBrick, robust=TRUE)
## Plotting the results
autoplot(sBrick)
autoplot(stlfBrick)

#The stlf function actually ends up creating less noise in the forecast and actually can be a better predictor of brick production.

#q7
## Plotting the dataset to evaluate which method to choose
autoplot(writing)
#i ended up choosing the drift method 
stlf(writing, method='rwdrift')

## Box Cox lambda
writingBC <- stlf(writing, method='rwdrift', robust=TRUE, lambda = BoxCox.lambda(writing))
autoplot(writingBC)


#8
autoplot(fancy)

#the fancy dataset shows a similar trend and seasonality but in a positive trend.

stlf(fancy,method='rwdrift')

fancyBC <- stlf(fancy, method='rwdrift', robust=TRUE, lambda = BoxCox.lambda(fancy))
autoplot(fancyBC)