Hyndman Excercises Ch. 5/6

Demetri Chokshi-Fox

Chapter 5

Question 1

library(fpp2)

## Warning: package 'fpp2' was built under R version 3.5.3

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.5.3

## Loading required package: forecast

## Warning: package 'forecast' was built under R version 3.5.3

## Loading required package: fma

## Warning: package 'fma' was built under R version 3.5.3

## Loading required package: expsmooth

## Warning: package 'expsmooth' was built under R version 3.5.3

daily20 = head(elecdaily, 20)
autoplot(daily20)

reg1 = tslm(Demand~ Temperature, data=daily20)
reg1

## 
## Call:
## tslm(formula = Demand ~ Temperature, data = daily20)
## 
## Coefficients:
## (Intercept)  Temperature  
##      39.212        6.757

1A) There is a positive relationship because as temperature rises, more people use air conditioning systems, driving up demand for electricity.

checkresiduals(reg1)

## 
##  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

1B) The model is adequate. The histogram is slightly skewed left, which may affect the prediction intervals. The autocorrelation plot shows a significant spike at lag 4 but it isn’t enough to the Breusch-Godfrey to be significant. Neither of these factors should affect the model’s viability.

newdata= data.frame(Temperature = 15)
fcast15 = forecast(reg1, h=1, newdata=newdata)
autoplot(fcast15)

newdata= data.frame(Temperature = 35)
fcast35 = forecast(reg1, h=1, newdata=newdata)
autoplot(fcast35)

1C) I believe that the forecast for 15 degrees is accurate but I don’t think that the forecast for 35 degrees is accurate. 35 degrees C is 95 degrees F, and that’s pretty toasty. I think that demand would be higher than forecasted.

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)))

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 3.857143       140.5701 108.6810 172.4591  90.21166 190.9285
## 4.000000       275.7146 245.2278 306.2014 227.57056 323.8586

sigma(fit)

## [1] 22.00373

#1D)
  #15 degrees:
#Upper Limit:
140.5701+(sigma(fit)*1.96)

## [1] 183.6974

#Lower Limit:
140.5701-(sigma(fit)*1.96)

## [1] 97.44279

#35 degrees:
#Upper Limit:
275.7146+(sigma(fit)*1.96)

## [1] 318.8419

#Lower limit:
  275.7146-(sigma(fit)*1.96)

## [1] 232.5873

1E)

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

This says my model is fairly accurate on the higher end of temperature, but less accurate on the lower end of temperature. This is most likely because the values in my sample were all higher temperatures, so the model was trained on that.

Question 2

autoplot(mens400)

2A)The function has an overall declining slope. There are also gaps in the function where there is missing data.

2B)

time.mens400 = time(mens400)
tslm.mens400 = tslm(mens400 ~ time.mens400, data= mens400)
autoplot(mens400)+ 
geom_abline(slope=tslm.mens400$coefficients[2], intercept=tslm.mens400$coefficients[1], color="red")

#Average rate of decrease in times per year:
  tslm.mens400$coefficents[2]

## NULL

cbind(Time = time.mens400, 
      Residuals = tslm.mens400$residuals) %>%
  as.data.frame() %>%
  ggplot(aes(x = Time, y = Residuals)) +
    geom_point() +
    ylab("Residuals of Regression Line")

## Warning: Removed 3 rows containing missing values (geom_point).

checkresiduals(tslm.mens400)

## 
##  Breusch-Godfrey test for serial correlation of order up to 6
## 
## data:  Residuals from Linear regression model
## LM test = 3.6082, df = 6, p-value = 0.7295

This says the fitted line is suitable, with one exception, the outlier in 1896. The checkresiduals function also indicates the same.

fcast.mens400 = forecast(tslm.mens400, newdata = data.frame(time.mens400 = 2020))
autoplot(mens400) + autolayer(fcast.mens400)

#Prediction Intervals
fcast.mens400$upper

## Time Series:
## Start = 2020 
## End = 2020 
## Frequency = 0.25 
##      Series 1 Series 2
## 2020 43.63487 44.53176

fcast.mens400$lower

## Time Series:
## Start = 2020 
## End = 2020 
## Frequency = 0.25 
##      Series 1 Series 2
## 2020 40.44975 39.55286

Question 3

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

The results are the counts of occurences of easter in each quarter from 1956 to 2010 for a quarterly australian beer sales dataset. Can be fractional if easter spans days that are in multiple quarters.

Question 5

autoplot(fancy)

The pattern is highly seasonal, with a large spike towards the end of the year and a small spike around March, corresponding with Christmas tourists and the local surfing festival, respectively. Each spike grows larger as time goes on.

5B Transforming the data into logs is necessary because it is clearly a nonlinear trend.

Time = time(fancy)
surffest = c()
for(i in 1:length(Time)){
  month = round(12*(Time[i]-floor(Time[i]))) + 1
  year = floor(Time[i])
  if(year>=1998 & month == 3){
  surffest[i] = 1 
  } else {
  surffest[i] = 0}
}
fancy.tslm = tslm(log(fancy)~ trend+season+surffest, data=fancy)

autoplot(fancy.tslm$residuals)

cbind(Residuals = fancy.tslm$residuals, 
  FittedValues = fancy.tslm$fitted.values)%>%
  as.data.frame%>%
  ggplot(aes(x=FittedValues, y=Residuals))+geom_point()

The residuals have a pattern with time, which means that the residuals have a correlation with time.

  cbind.data.frame(
    Month = factor(
      month.abb[round(12*(Time - floor(Time)) + 1)],
      labels = month.abb,
      ordered = TRUE
    ),
    Residuals = fancy.tslm$residuals
    ) %>%
  ggplot(aes(x = Month,
             y = Residuals)) +
    geom_boxplot()

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

The box plots show a large range of residuals in August, but no large problems.

fancy.tslm$coefficients

## (Intercept)       trend     season2     season3     season4     season5 
##  7.60586042  0.02239298  0.25104367  0.69520658  0.38293405  0.40799437 
##     season6     season7     season8     season9    season10    season11 
##  0.44696253  0.60821562  0.58535238  0.66634464  0.74403359  1.20301639 
##    season12    surffest 
##  1.95813657          NA

The coefficients tell me that the largest effects are in the later months, as well as in March (season3). There is a distictly increasing trend, however it is not a large yearly increase (2%).

#fancy.residuals <- checkresiduals(fancy.tslm)
#fancy.residuals
#returns error:
 # Error in chol2inv(auxfit$qr$qr) : 
  #element (31, 31) is zero, so the inverse cannot be computed

There are significant spikes at lags 1, 2, 10, 21, 23, 24, 26, and 28, indicating autocorrelation.

newdata= data.frame(sample = rep(0,84))
threeYrFcast = as.data.frame(forecast(fancy.tslm, h=36, newdata=newdata))

## Warning in forecast.lm(fancy.tslm, h = 36, newdata = newdata): Could not
## find required variable surffest in newdata. Specify newdata as a named
## data.frame

## Warning in predict.lm(predict_object, newdata = newdata, se.fit = TRUE, :
## prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(predict_object, newdata = newdata, se.fit = TRUE, :
## prediction from a rank-deficient fit may be misleading

threeYrFcast

##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## Jan 1994       9.509264  9.246995  9.771532  9.105002  9.913525
## Feb 1994       9.782700  9.520432 10.044969  9.378439 10.186962
## Mar 1994      10.249256  9.986988 10.511525  9.844995 10.653518
## Apr 1994       9.959377  9.697108 10.221645  9.555115 10.363638
## May 1994      10.006830  9.744562 10.269098  9.602569 10.411091
## Jun 1994      10.068191  9.805923 10.330459  9.663930 10.472452
## Jul 1994      10.251837  9.989569 10.514105  9.847576 10.656099
## Aug 1994      10.251367  9.989099 10.513635  9.847106 10.655628
## Sep 1994      10.354752 10.092484 10.617020  9.950491 10.759014
## Oct 1994      10.454834 10.192566 10.717102 10.050573 10.859095
## Nov 1994      10.936210 10.673942 11.198478 10.531949 11.340471
## Dec 1994      11.713723 11.451455 11.975991 11.309462 12.117984
## Jan 1995       9.777979  9.512777 10.043182  9.369195 10.186763
## Feb 1995      10.051416  9.786214 10.316618  9.642632 10.460200
## Mar 1995      10.517972 10.252770 10.783174 10.109188 10.926756
## Apr 1995      10.228092  9.962890 10.493295  9.819308 10.636876
## May 1995      10.275546 10.010343 10.540748  9.866762 10.684330
## Jun 1995      10.336907 10.071704 10.602109  9.928123 10.745691
## Jul 1995      10.520553 10.255351 10.785755 10.111769 10.929337
## Aug 1995      10.520083 10.254880 10.785285 10.111299 10.928867
## Sep 1995      10.623468 10.358266 10.888670 10.214684 11.032252
## Oct 1995      10.723550 10.458347 10.988752 10.314766 11.132334
## Nov 1995      11.204926 10.939723 11.470128 10.796142 11.613710
## Dec 1995      11.982439 11.717236 12.247641 11.573655 12.391223
## Jan 1996      10.046695  9.777950 10.315440  9.632451 10.460940
## Feb 1996      10.320132 10.051387 10.588877  9.905887 10.734376
## Mar 1996      10.786688 10.517943 11.055433 10.372443 11.200932
## Apr 1996      10.496808 10.228063 10.765553 10.082564 10.911053
## May 1996      10.544261 10.275516 10.813007 10.130017 10.958506
## Jun 1996      10.605623 10.336878 10.874368 10.191378 11.019867
## Jul 1996      10.789269 10.520524 11.058014 10.375024 11.203513
## Aug 1996      10.788798 10.520053 11.057543 10.374554 11.203043
## Sep 1996      10.892184 10.623439 11.160929 10.477939 11.306428
## Oct 1996      10.992266 10.723521 11.261011 10.578021 11.406510
## Nov 1996      11.473641 11.204896 11.742386 11.059397 11.887886
## Dec 1996      12.251155 11.982409 12.519900 11.836910 12.665399
## Jan 1997      10.315411 10.042538 10.588283  9.894804 10.736018
## Feb 1997      10.588848 10.315975 10.861720 10.168241 11.009454
## Mar 1997      11.055403 10.782531 11.328276 10.634797 11.476010
## Apr 1997      10.765524 10.492651 11.038396 10.344917 11.186131
## May 1997      10.812977 10.540105 11.085850 10.392371 11.233584
## Jun 1997      10.874338 10.601466 11.147211 10.453732 11.294945
## Jul 1997      11.057984 10.785112 11.330857 10.637378 11.478591
## Aug 1997      11.057514 10.784642 11.330387 10.636907 11.478121
## Sep 1997      11.160899 10.888027 11.433772 10.740293 11.581506
## Oct 1997      11.260981 10.988109 11.533854 10.840375 11.681588
## Nov 1997      11.742357 11.469485 12.015230 11.321750 12.162964
## Dec 1997      12.519870 12.246998 12.792743 12.099264 12.940477
## Jan 1998      10.584127 10.306568 10.861685 10.156297 11.011957
## Feb 1998      10.857563 10.580005 11.135122 10.429733 11.285393
## Mar 1998      11.324119 11.046561 11.601678 10.896289 11.751949
## Apr 1998      11.034240 10.756681 11.311798 10.606410 11.462070
## May 1998      11.081693 10.804134 11.359252 10.653863 11.509523
## Jun 1998      11.143054 10.865495 11.420613 10.715224 11.570884
## Jul 1998      11.326700 11.049142 11.604259 10.898870 11.754530
## Aug 1998      11.326230 11.048671 11.603789 10.898400 11.754060
## Sep 1998      11.429615 11.152056 11.707174 11.001785 11.857445
## Oct 1998      11.529697 11.252138 11.807256 11.101867 11.957527
## Nov 1998      12.011073 11.733514 12.288632 11.583243 12.438903
## Dec 1998      12.788586 12.511027 13.066145 12.360756 13.216416
## Jan 1999      10.852842 10.570067 11.135618 10.416971 11.288714
## Feb 1999      11.126279 10.843503 11.409055 10.690407 11.562151
## Mar 1999      11.592835 11.310059 11.875611 11.156963 12.028707
## Apr 1999      11.302955 11.020180 11.585731 10.867084 11.738827
## May 1999      11.350409 11.067633 11.633185 10.914537 11.786280
## Jun 1999      11.411770 11.128994 11.694546 10.975898 11.847642
## Jul 1999      11.595416 11.312640 11.878192 11.159544 12.031288
## Aug 1999      11.594946 11.312170 11.877722 11.159074 12.030817
## Sep 1999      11.698331 11.415555 11.981107 11.262459 12.134203
## Oct 1999      11.798413 11.515637 12.081189 11.362541 12.234285
## Nov 1999      12.279789 11.997013 12.562564 11.843917 12.715660
## Dec 1999      13.057302 12.774526 13.340078 12.621430 13.493174
## Jan 2000      11.121558 10.833063 11.410053 10.676871 11.566246
## Feb 2000      11.394995 11.106500 11.683490 10.950307 11.839682
## Mar 2000      11.861551 11.573056 12.150046 11.416863 12.306238
## Apr 2000      11.571671 11.283176 11.860166 11.126984 12.016359
## May 2000      11.619124 11.330629 11.907620 11.174437 12.063812
## Jun 2000      11.680486 11.391991 11.968981 11.235798 12.125173
## Jul 2000      11.864132 11.575637 12.152627 11.419444 12.308819
## Aug 2000      11.863661 11.575166 12.152157 11.418974 12.308349
## Sep 2000      11.967047 11.678552 12.255542 11.522359 12.411734
## Oct 2000      12.067129 11.778633 12.355624 11.622441 12.511816
## Nov 2000      12.548504 12.260009 12.837000 12.103817 12.993192
## Dec 2000      13.326018 13.037522 13.614513 12.881330 13.770705

#intervals
upper = threeYrFcast$upper
lower = threeYrFcast$lower
predictionIntervals = as.data.frame(upper, lower)
predictionIntervals

## data frame with 0 columns and 0 rows

pred.intvl.invlog = exp(predictionIntervals)
pred.intvl.invlog

## data frame with 0 columns and 0 rows

I could probably improve these predictions by making the data a log-log function.

Question 6

str(gasoline)

##  Time-Series [1:1355] from 1991 to 2017: 6.62 6.43 6.58 7.22 6.88 ...

gasoline2004 = window(gasoline, end=2004)
reg6.1 = tslm(gasoline2004~trend + fourier(gasoline2004, K=2))
reg6.2 = tslm(gasoline2004~trend + fourier(gasoline2004, K=5))
reg6.3 = tslm(gasoline2004~trend + fourier(gasoline2004, K=7))
reg6.4 = tslm(gasoline2004~trend + fourier(gasoline2004, K=10))
plot.ts(gasoline2004, reg6.1$fitted.values)

plot.ts(gasoline2004, reg6.2$fitted.values)

plot.ts(gasoline2004, reg6.3$fitted.values)

plot.ts(gasoline2004, reg6.4$fitted.values)

The plots all show high correlation.

AIC(reg6.1)

## [1] 247.8825

AIC(reg6.2)

## [1] 202.2442

AIC(reg6.3)

## [1] 166.4069

AIC(reg6.4)

## [1] 167.2706

reg6.3 with k=7 is the best fit.

checkresiduals(reg6.3)

## 
##  Breusch-Godfrey test for serial correlation of order up to 104
## 
## data:  Residuals from Linear regression model
## LM test = 145.31, df = 104, p-value = 0.004671

fc = forecast(reg6.3, newdata=data.frame(fourier(gasoline2004,K=7,h=52)))
plot(fc)

autoplot(fc) +autolayer(window(gasoline2004))

Question 7

autoplot(huron)

h <- 10
reg7 = tslm(huron~trend, data=huron)
time = time(huron)
time.break = 1915
tb1 = ts(pmax(0, time - time.break), start = 1875)
fit.pw = tslm(huron~time+tb1)
t.new = time[length(time)] + seq(h)
tb1.new = tb1[length(tb1)] + seq(h)
newdata = cbind(time=t.new, tb1=tb1.new)%>%
as.data.frame()

reg7Fcast = forecast(reg7)
reg7Fcast

##      Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## 1973       7.806127 6.317648 9.294605 5.516501 10.095752
## 1974       7.781926 6.292536 9.271315 5.490899 10.072952
## 1975       7.757724 6.267406 9.248043 5.465269 10.050179
## 1976       7.733523 6.242259 9.224788 5.439613 10.027434
## 1977       7.709322 6.217094 9.201550 5.413929 10.004715
## 1978       7.685121 6.191912 9.178331 5.388219  9.982024
## 1979       7.660920 6.166712 9.155128 5.362481  9.959359
## 1980       7.636719 6.141494 9.131943 5.336717  9.936721
## 1981       7.612518 6.116259 9.108776 5.310925  9.914110
## 1982       7.588317 6.091007 9.085626 5.285107  9.891526

reg7.1Fcast = forecast(fit.pw, newdata = newdata)
reg7.1Fcast

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1973       8.455119 7.063583 9.846655 6.314483 10.59575
## 1974       8.454992 7.061518 9.848467 6.311374 10.59861
## 1975       8.454866 7.059398 9.850333 6.308182 10.60155
## 1976       8.454739 7.057225 9.852253 6.304906 10.60457
## 1977       8.454612 7.054998 9.854227 6.301549 10.60768
## 1978       8.454486 7.052717 9.856254 6.298109 10.61086
## 1979       8.454359 7.050384 9.858334 6.294587 10.61413
## 1980       8.454232 7.047997 9.860467 6.290984 10.61748
## 1981       8.454106 7.045558 9.862653 6.287300 10.62091
## 1982       8.453979 7.043067 9.864891 6.283535 10.62442

autoplot(huron)+
autolayer(fitted(reg7), series = "Linear") +
  autolayer(fitted(fit.pw), series = "Piecewise") +
  autolayer(reg7.1Fcast, series="Piecewise") +
  autolayer(reg7Fcast, series="Linear", PI=FALSE) +
  xlab("Year") + ylab("Water Level") +
  ggtitle("Water level Change") +
  guides(colour = guide_legend(title = " "))

. . .

Chapter 6

library(fpp2)
plot(plastics)

There is an increasing trend with obvious seasonality.

multdecomp = decompose(plastics, type = "multiplicative")
plot(multdecomp)

2C Yes, the results show a clear increasing trend, seasonality, and cyclical behavior.

library(seasonal)

## Warning: package 'seasonal' was built under R version 3.5.3

autoplot(plastics, series="Data") +
  autolayer(trendcycle(multdecomp), series="Trend") +
  autolayer(seasadj(multdecomp), series="Seasonally Adjusted") +
  xlab("Year") + ylab("Sales in $") +
  ggtitle("Sales of Product A") +
  scale_colour_manual(values=c("gray","blue","red"),
                      breaks=c("Data","Seasonally Adjusted","Trend"))

## Warning: Removed 12 rows containing missing values (geom_path).

plastics[15] = plastics[15] + 500
multdecomp1 = decompose(plastics, type= "multiplicative")
plot(multdecomp1)

2E. The outlier adds a blip to the seasonality plot and also increases the trend at a faster rate through the first three periods on the trend plot.

2F. Yes, it makes a difference in the trend line.

Can’t do Q4, don;t have the data

Question 4

There is a clear upwards trend. There is clear seasonality, although it doesn’t have a clear meaning when plotted monthly. The remainder in 1991 is unexpectedly large, indicating a possible recession.

Question 5

autoplot(cangas)

ggsubseriesplot(cangas)

ggseasonplot(cangas)

The population growth drives the growth. Seasonality is driven by seasonal demand, with demand spiking in extreme temperature ranges.

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

cangasx11 = seas(cangas, x11="")
autoplot(cangasx11)

cangas.seas = seas(cangas)
autoplot(cangas.seas)

The X11 decomp shows a smoother trend and less fluctuation, which I interpret as a better decomposition, as it accounts for more of the data.

Question 6

stlfixed = stl(bricksq, s.window = "periodic", robust = TRUE)
autoplot(stlfixed)

stlchanging = stl(bricksq, s.window = 13)
autoplot(stlchanging)

fixed.bricksq.fc = forecast(stlchanging)
fixed.bricksq.fc

##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1994 Q4       465.7326 437.4903 493.9748 422.5398 508.9254
## 1995 Q1       424.6727 384.7119 464.6335 363.5580 485.7874
## 1995 Q2       483.9914 435.0232 532.9595 409.1010 558.8817
## 1995 Q3       493.9988 437.4242 550.5734 407.4754 580.5222
## 1995 Q4       465.7326 402.4454 529.0198 368.9431 562.5220
## 1996 Q1       424.6727 355.3067 494.0387 318.5865 530.7589
## 1996 Q2       483.9914 409.0259 558.9568 369.3416 598.6411
## 1996 Q3       493.9988 413.8128 574.1848 371.3650 616.6326

autoplot(fixed.bricksq.fc)

stlfixed %>% seasadj() %>% naive() %>%
  autoplot() + ylab("New orders index") +
  ggtitle("Naive forecasts of seasonally adjusted data")

bricksq.fs = stlf(bricksq)
autoplot(bricksq.fs)

checkresiduals(bricksq.fs)

## Warning in checkresiduals(bricksq.fs): The fitted degrees of freedom is
## based on the model used for the seasonally adjusted data.

## 
##  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

brickq.robust = stlf(bricksq, robust = TRUE)
autoplot(brickq.robust)

checkresiduals(brickq.robust)

## Warning in checkresiduals(brickq.robust): The fitted degrees of freedom is
## based on the model used for the seasonally adjusted data.

## 
##  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

b.train <- subset(bricksq,end = length(bricksq) - 8)
b.test <- subset(bricksq,start = length(bricksq) - 7)
b.train.snaive <- snaive(b.train)
fs1 <- forecast(b.train.snaive, data=b.test)
plot(fs1)

b.train.stlf = stlf(b.train, robust = TRUE)
fs2 = forecast(b.train.stlf, data= b.test)
plot(fs2)

Question 7

autoplot(writing)

fcast = stlf(writing, method = 'naive')
autoplot(fcast)

fcast1 = stlf(writing, method = 'rwdrift')
autoplot(fcast1)

fcast2 = stlf(writing, method = 'naive', lambda = BoxCox.lambda(writing))
autoplot(fcast2)

fcast3 = stlf(writing, method='rwdrift', lambda = BoxCox.lambda(writing))
autoplot(fcast3)

Question 8

``` autoplot(fancy) fcast.fancy = stlf(fancy, method=‘naive’) autoplot(fcast.fancy)

fcast1.fancy = stlf(fancy, method=‘rwdrift’) autoplot(fcast1.fancy)

fcast2.fancy = stlf(fancy, method=‘naive’, lambda = BoxCox.lambda(fancy)) autoplot(fcast2.fancy)

fcast3.fancy = stlf(fancy, method=‘rwdrift’, lambda = BoxCox.lambda(fancy)) autoplot(fcast3.fancy)