Assignment1

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

#Chapter5 Exercise

#5.1 Daily Electricty Demand 
daily20<- head(elecdaily,20)
#a
plot(daily20)

We notice that there is positive relationship between demand and temperature since the more air conditioner would be used when the temperature goes high, especially in summer.

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

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

#c
summary(fit)

## 
## Call:
## tslm(formula = Demand ~ Temperature, data = daily20)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.060  -7.117  -1.437  17.484  27.102 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.2117    17.9915   2.179   0.0428 *  
## Temperature   6.7572     0.6114  11.052 1.88e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22 on 18 degrees of freedom
## Multiple R-squared:  0.8716, Adjusted R-squared:  0.8644 
## F-statistic: 122.1 on 1 and 18 DF,  p-value: 1.876e-09

Min Tem= 15, demand= 142 Max Tem= 35, demand= 276

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

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

The overall dataset for ‘elecdaily’ shows that the linear model which is estimated by first 20 entries is not precise. The demand of electricity would decrease around 20. There’re two situations that the usage of electricity would goes up when weather is hot or cold.

#5.2 Olympic Games
#a 
autoplot(mens400)

The winning times of men’s 400 meters constantly decreased in years.

autoplot(mens400)+
  geom_point() +
  geom_smooth(method="lm", se=FALSE)

## Warning: Removed 3 rows containing non-finite values (stat_smooth).

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

#c
fit2<-tslm(mens400~trend)
print(fit2)

## 
## Call:
## tslm(formula = mens400 ~ trend)
## 
## Coefficients:
## (Intercept)        trend  
##     50.3078      -0.2583

checkresiduals(fit2)

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

#d
rwf(mens400, drift=TRUE, h=1)

##      Point Forecast   Lo 80    Hi 80    Lo 95   Hi 95
## 2020        42.6163 41.0123 44.22029 40.16319 45.0694

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

This dataset describe the exact date of the Easter holiday. Usually, Easter holiday would belongs to March(Q1=1) or April(Q2=1), sometimes it seperates between April and March like 1/3, and 2/3.

#5.5
#a
plot(fancy)

There is seasonal fluctuation, the sales would goes up at the end of years.

#b
fancy1<-as.data.frame(fancy)


fancy2<-as.numeric(unlist(fancy1))

p<-hist(fancy2,col='blue')

fancy1$logx<-log(fancy1$x)

The log transformation of dataset would normalize the distribution

#c

fit5<-tslm(BoxCox(fancy,0)~trend+season)
summary(fit5)

## 
## Call:
## tslm(formula = BoxCox(fancy, 0) ~ trend + season)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41644 -0.12619  0.00608  0.11389  0.38567 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.6058604  0.0768740  98.939  < 2e-16 ***
## trend       0.0223930  0.0008448  26.508  < 2e-16 ***
## season2     0.2510437  0.0993278   2.527 0.013718 *  
## season3     0.6952066  0.0993386   6.998 1.18e-09 ***
## season4     0.3829341  0.0993565   3.854 0.000252 ***
## season5     0.4079944  0.0993817   4.105 0.000106 ***
## season6     0.4469625  0.0994140   4.496 2.63e-05 ***
## season7     0.6082156  0.0994534   6.116 4.69e-08 ***
## season8     0.5853524  0.0995001   5.883 1.21e-07 ***
## season9     0.6663446  0.0995538   6.693 4.27e-09 ***
## season10    0.7440336  0.0996148   7.469 1.61e-10 ***
## season11    1.2030164  0.0996828  12.068  < 2e-16 ***
## season12    1.9581366  0.0997579  19.629  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1858 on 71 degrees of freedom
## Multiple R-squared:  0.9527, Adjusted R-squared:  0.9447 
## F-statistic: 119.1 on 12 and 71 DF,  p-value: < 2.2e-16

#d
checkresiduals(fit5)

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

#e
box <-residuals(fit5)

box1<-data.frame(box,Month=rep(1:12,7))
boxplot(box~Month,data=box1)

The sales experienced most volatile fluctuations on March.

#g
bg<-bgtest(fit5, order=1, data=fancy2)
bg

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  fit5
## LM test = 21.54, df = 1, p-value = 3.466e-06

#h
rwf(fancy, drift=TRUE, h=36)

##          Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
## Jan 1994       105901.6 87650.115 124153.1  77988.371 133814.8
## Feb 1994       107142.5 81025.551 133259.4  67200.075 147084.9
## Mar 1994       108383.4 76027.016 140739.8  58898.574 157868.2
## Apr 1994       109624.3 71840.195 147408.5  51838.485 167410.2
## May 1994       110865.2 68154.500 153576.0  45544.800 176185.7
## Jun 1994       112106.2 64813.059 159399.2  39777.607 184434.7
## Jul 1994       113347.1 61724.118 164970.0  34396.579 192297.6
## Aug 1994       114588.0 58828.823 170347.1  29311.707 199864.3
## Sep 1994       115828.9 56086.940 175570.9  24461.459 207196.3
## Oct 1994       117069.8 53469.641 180670.0  19801.746 214337.9
## Nov 1994       118310.7 50955.514 185665.9  15299.820 221321.6
## Dec 1994       119551.6 48528.192 190575.1  10930.651 228172.6
## Jan 1995       120792.6 46174.873 195410.2   6674.660 234910.4
## Feb 1995       122033.5 43885.344 200181.6   2516.228 241550.7
## Mar 1995       123274.4 41651.329 204897.4  -1557.304 248106.1
## Apr 1995       124515.3 39466.022 209564.6  -5556.343 254586.9
## May 1995       125756.2 37323.762 214188.7  -9489.547 261002.0
## Jun 1995       126997.1 35219.788 218774.5 -13364.198 267358.4
## Jul 1995       128238.0 33150.059 223326.0 -17186.475 273662.5
## Aug 1995       129478.9 31111.119 227846.8 -20961.665 279919.6
## Sep 1995       130719.9 29099.985 232339.7 -24694.329 286134.1
## Oct 1995       131960.8 27114.071 236807.5 -28388.423 292310.0
## Nov 1995       133201.7 25151.116 241252.3 -32047.404 298450.8
## Dec 1995       134442.6 23209.137 245676.1 -35674.304 304559.5
## Jan 1996       135683.5 21286.383 250080.7 -39271.802 310638.8
## Feb 1996       136924.4 19381.302 254467.6 -42842.272 316691.1
## Mar 1996       138165.3 17492.512 258838.2 -46387.829 322718.5
## Apr 1996       139406.3 15618.776 263193.7 -49910.360 328722.9
## May 1996       140647.2 13758.987 267535.4 -53411.563 334705.9
## Jun 1996       141888.1 11912.145 271864.0 -56892.964 340669.1
## Jul 1996       143129.0 10077.350 276180.7 -60355.940 346613.9
## Aug 1996       144369.9  8253.785 280486.0 -63801.742 352541.6
## Sep 1996       145610.8  6440.709 284781.0 -67231.504 358453.2
## Oct 1996       146851.7  4637.445 289066.0 -70646.258 364349.7
## Nov 1996       148092.7  2843.376 293341.9 -74046.950 370232.3
## Dec 1996       149333.6  1057.938 297609.2 -77434.441 376101.6

#5.6 Gasoline
#a 
weeklygas<-window(gasoline, end=2005)

length(weeklygas)

## [1] 726

weeklygas1<- tslm(weeklygas~ trend + fourier(weeklygas, K=26))

autoplot(weeklygas, series= 'Observed Gasoline')+autolayer(fitted(weeklygas1), series='Fitted Value') +
ggtitle("K=26")+xlab('Year')+ylab ('Gasoline')

#K=9
weeklygas2<- tslm(weeklygas~ trend + fourier(weeklygas, K=9))

autoplot(weeklygas, series= 'Observed Gasoline')+autolayer(fitted(weeklygas2), series='Fitted Value') +
  ggtitle("K=9")+xlab('Year')+ylab ('Gasoline')

#K=4
weeklygas3<- tslm(weeklygas~ trend + fourier(weeklygas, K=4))

autoplot(weeklygas, series= 'Observed Gasoline')+autolayer(fitted(weeklygas3), series='Fitted Value') +
  ggtitle("K=4")+xlab('Year')+ylab ('Gasoline')

#b
CV(weeklygas1)

##            CV           AIC          AICc           BIC         AdjR2 
##  7.426757e-02 -1.889694e+03 -1.880500e+03 -1.637379e+03  8.526098e-01

CV(weeklygas2)

##            CV           AIC          AICc           BIC         AdjR2 
##  7.161451e-02 -1.912292e+03 -1.910980e+03 -1.815954e+03  8.506543e-01

CV(weeklygas3)

##            CV           AIC          AICc           BIC         AdjR2 
##  7.648608e-02 -1.864112e+03 -1.863742e+03 -1.813649e+03  8.382404e-01

We compared three harmonic regression models, and select one with lowest AIC&CV, and highest adjusted R^2: K=9

#c
checkresiduals(weeklygas2)

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

#d

wg<-forecast(weeklygas2, newdata=data.frame(fourier(weeklygas,9,52)))

weeklygasoline<-window(gasoline, start=2005, end=2006)

autoplot(weeklygasoline, series=' Observed Gasoline') + autolayer(wg$mean, series='Fitted Value')+
  ggtitle("Comparing Prediction and Actual")+xlab('Year')+ylab('Gasoline')

#5.7
#a
autoplot(huron)+ xlab('Year')+ ggtitle('Water Level of Lake Huron')

The water level of lake has seasonal trend, there is severe fluctuations periodically.

#b
#Linear Reg
fit7<-tslm(huron~trend)
summary(fit7)

## 
## Call:
## tslm(formula = huron ~ trend)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50997 -0.72726  0.00083  0.74402  2.53565 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.202037   0.230111  44.335  < 2e-16 ***
## trend       -0.024201   0.004036  -5.996 3.55e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.13 on 96 degrees of freedom
## Multiple R-squared:  0.2725, Adjusted R-squared:  0.2649 
## F-statistic: 35.95 on 1 and 96 DF,  p-value: 3.545e-08

#Piecewise Linear Trend
t <- time(huron)
t1 <- ts(pmax(0, t - 1915), start = 1875)
fitpw <- tslm(huron ~ t + t1)
summary(fitpw)

## 
## Call:
## tslm(formula = huron ~ t + t1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.49626 -0.66240 -0.07139  0.85163  2.39222 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 132.90870   19.97687   6.653 1.82e-09 ***
## t            -0.06498    0.01051  -6.181 1.58e-08 ***
## t1            0.06486    0.01563   4.150 7.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.045 on 95 degrees of freedom
## Multiple R-squared:  0.3841, Adjusted R-squared:  0.3711 
## F-statistic: 29.62 on 2 and 95 DF,  p-value: 1.004e-10

#c
forecast(fit7, h=8)

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

#Chapter 6 Exercise
#6.1 Plastic Manufacturers
#2a
autoplot(plastics)+ ggtitle('The Monthly Sales(in thousands) of Product A')

There is an obvious seasonality, the sales would reachs the annual peak in Q3. And overall graph exhibits positive trend.

#2b
plastics %>% decompose(type="multiplicative") %>% 
  autoplot() + xlab("Year") +
  ggtitle("Classical multiplicative decomposition 
          of plastics")

#2c The result of multiplicaetive decomposition is consistent with the graph in part a.

#2d

multi<-decompose(plastics, type="multiplicative")

seasonal_adj<-multi$x/multi$seasonal
autoplot(seasonal_adj)+ ylab('Seasonally Adjusted Data')

#2e
pl<-plastics
pl[32]<-pl[32]+500

multi1<-decompose(pl, type="multiplicative")
seasonal_adj1<-multi1$x/multi1$seasonal
autoplot(seasonal_adj1)+ ggtitle('Seasonally Adjusted Data(Outlier Near The Middle')

#2f
pl1<-plastics
pl1[57]<-pl[57]+500

multi2<-decompose(pl1, type="multiplicative")
seasonal_adj2<-multi2$x/multi1$seasonal
autoplot(seasonal_adj2)+ ggtitle('Seasonally Adjusted Data(Outlier Near The End')

#6.3 Decomposition
retaildata<- readxl::read_xlsx("retail.xlsx",skip=1)

myts <- ts(retaildata[,"A3349873A"],
           frequency=12, start=c(1982,4))

myts %>% seas(x11="") %>%
  autoplot() +
  ggtitle("X11 decomposition of Retail Data")

6.4 Monthly civilian labour force in Australia a: It has increasing trends in seasonal adjustment plot. In the remainder part, some months between 1991 to 1992 experienced flutuactions. b: Yes, the reduction between 1991 to 1992 up to 400 units.

#6.5 Monthly Canadian Gas Production
#5a
autoplot(cangas)

ggsubseriesplot(cangas)

ggseasonplot(cangas)

#5b
cangas %>%
  stl(t.window=46, s.window="periodic", robust=TRUE) %>%
  autoplot()+ ggtitle('STL decomposition of Monthly Canadian Gas Production 1.0 ')

cangas %>%
  stl(t.window=46, s.window=99, robust=TRUE) %>%
  autoplot()+ ggtitle('STL decomposition of Monthly Canadian Gas Production 2.0')

cangas %>%
  stl(t.window=46, s.window=30, robust=TRUE) %>%
  autoplot()+ ggtitle('STL decomposition of Monthly Canadian Gas Production 3.0')

cangas %>%
  stl(t.window=46, s.window=7, robust=TRUE) %>%
  autoplot()+ ggtitle('STL decomposition of Monthly Canadian Gas Production 4.0')

#5c
cangas %>% seas() %>%
  autoplot() +
  ggtitle("SEATS decomposition of Monthly Canadian Gas Production ")

cangas %>% seas(x11="") %>%
  autoplot() +
  ggtitle("X11 Monthly Canadian Gas Production")

#6.6 Australian Quarterly Clay Brick Production
#6a
bricksq %>%
  stl(t.window=39, s.window="periodic", robust=TRUE) %>%
  autoplot()+ 
  ggtitle('STL decomposition of Australian Quarterly Clay Brick Production(Fixed Seasonality)')

bricksq %>%
  stl(t.window=39, s.window=7, robust=TRUE) %>%
  autoplot()+ 
  ggtitle('STL decomposition of Australian Quarterly Clay Brick Production(Changing Seasonality)')

#6b
bricksq %>%
  stl(t.window=39, s.window="periodic", robust=FALSE) %>% seasadj() %>%
  autoplot()+
  ggtitle('STL decomposition of Australian Quarterly Clay Brick Production(Seasonally Adjusted)')

#6c
bsq <- stl(bricksq, t.window=39, s.window="periodic", robust=FALSE)

bsq %>% seasadj() %>% naive() %>% autoplot() +
  ggtitle("Naive forecasts of seasonally adjusted data")

#6d
rwf(stlf(bricksq,method='naive'),h=8)

##         Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## 1996 Q4       494.0000 394.8325 593.1675 342.33640 645.6636
## 1997 Q1       421.7450 281.5010 561.9891 207.26030 636.2298
## 1997 Q2       566.2550 394.4918 738.0181 303.56590 828.9440
## 1997 Q3       383.4956 185.1606 581.8306  80.16835 686.8228
## 1997 Q4       604.5044 382.7592 826.2497 265.37431 943.6346
## 1998 Q1       494.0000 251.0902 736.9098 122.50156 865.4984
## 1998 Q2       421.7450 159.3725 684.1176  20.48085 823.0092
## 1998 Q3       566.2550 285.7669 846.7430 137.28552 995.2244

#6e
checkresiduals(stlf(bricksq,method='naive'))

## Warning in checkresiduals(stlf(bricksq, method = "naive")): The fitted
## degrees of freedom is based on the model used for the seasonally adjusted
## data.

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  Random walk
## Q* = 40.829, df = 8, p-value = 2.244e-06
## 
## Model df: 0.   Total lags used: 8

#6f

bsq1 <- stl(bricksq, t.window=39, s.window="periodic", robust=TRUE)

bsq1 %>% seasadj() %>% naive() %>% autoplot() +
  ggtitle("Robust STL decomposition of seasonally adjusted data")+
  xlab('Year')+ ylab('Clay Brick Production')

#6g
bsqtrain<-window(bricksq, end=c(1991,4))
bsqtest<-as.numeric(window(bricksq,start= c(1992,1), end= c(1994,3)))

bsqna<-snaive(bsqtrain)

bsqnaf<-rwf(bsqna, h=11)
navr <- bsqtest- bsqnaf$mean


bsqstlf<-stlf(bsqtrain)
bsqstlff<-rwf(bsqstlf, h=11)
stlfr<-bsqtest-bsqstlff$mean


autoplot(navr,series='Naive Method') + autolayer(stlfr, series='STLF')+ 
      ggtitle('Compare Forecasts')+
      xlab('Year')

#6.7 Writing Series
plot(writing)

hist(writing)

wrna <- stlf(writing, method='naive')
wrrwf <- stlf(writing, method='rwdrift')

wrnar<-as.numeric(residuals(wrna))
wrrwfr<-as.numeric(residuals(wrrwf))

wrnar<-na.omit(wrnar)
wrrwfr<-na.omit(wrrwfr)

rmsena<-sqrt(sum((wrnar)^2)/120)
rmserwf<-sqrt(sum((wrrwfr)^2)/120)

rmsena

## [1] 46.49735

rmserwf

## [1] 46.36373

#6.8 Fancy Series
fancyts=ts(fancy,frequency=12)
plot(fancyts)

hist(fancyts)

ffcna <- stlf(writing, method='naive')
ffcdrf <- stlf(writing, method='rwdrift')
ffclog<- stlf(writing, method='rwdrift', lambda=0)

ffcna

##          Point Forecast    Lo 80     Hi 80     Lo 95     Hi 95
## Jan 1978       990.5776 930.7390 1050.4162 899.06240 1082.0929
## Feb 1978      1021.2781 936.6535 1105.9026 891.85600 1150.7002
## Mar 1978      1075.1129 971.4694 1178.7564 916.60387 1233.6219
## Apr 1978      1002.2253 882.5481 1121.9025 819.19485 1185.2558
## May 1978       967.7282 833.9251 1101.5314 763.09393 1172.3625
## Jun 1978      1032.4996 885.9256 1179.0736 808.33399 1256.6652
## Jul 1978       881.2474 722.9294 1039.5655 639.12088 1123.3740
## Aug 1978       469.7581 300.5090  639.0072 210.91393  728.6023
## Sep 1978       909.9409 730.4251 1089.4567 635.39518 1184.4866
## Oct 1978       995.4665 806.2402 1184.6927 706.06988 1284.8630
## Nov 1978       929.5874 731.1252 1128.0495 626.06566 1233.1090
## Dec 1978       993.7330 786.4460 1201.0200 676.71493 1310.7511
## Jan 1979       990.5776 774.8265 1206.3288 660.61476 1320.5405
## Feb 1979      1021.2781 797.3826 1245.1736 678.85943 1363.6967
## Mar 1979      1075.1129 843.3590 1306.8668 720.67593 1429.5499
## Apr 1979      1002.2253 762.8709 1241.5797 636.16438 1368.2863
## May 1979       967.7282 721.0074 1214.4491 590.40124 1345.0552
## Jun 1979      1032.4996 778.6260 1286.3733 644.23336 1420.7659
## Jul 1979       881.2474 620.4170 1142.0778 482.34177 1280.1531
## Aug 1979       469.7581 202.1518  737.3644  60.48953  879.0267
## Sep 1979       909.9409 635.7260 1184.1558 490.56540 1329.3164
## Oct 1979       995.4665 714.7986 1276.1344 566.22196 1424.7110
## Nov 1979       929.5874 642.6115 1216.5632 490.69570 1368.4790
## Dec 1979       993.7330 700.5849 1286.8811 545.40174 1442.0643

ffcdrf

##          Point Forecast    Lo 80     Hi 80     Lo 95     Hi 95
## Jan 1978       994.1148 933.9446 1054.2850 902.09247 1086.1371
## Feb 1978      1028.3524 942.5528 1114.1521 897.13317 1159.5717
## Mar 1978      1085.7244 979.7839 1191.6649 923.70237 1247.7465
## Apr 1978      1016.3740 893.0618 1139.6862 827.78432 1204.9637
## May 1978       985.4141 846.4570 1124.3711 772.89758 1197.9306
## Jun 1978      1053.7226 900.3182 1207.1271 819.11078 1288.3345
## Jul 1978       906.0076 739.0422 1072.9731 650.65599 1161.3592
## Aug 1978       498.0555 318.2147  677.8962 223.01282  773.0981
## Sep 1978       941.7754 749.6073 1133.9435 647.87963 1235.6712
## Oct 1978      1030.8382 826.7912 1234.8851 718.77523 1342.9011
## Nov 1978       968.4962 752.9447 1184.0477 638.83869 1298.1537
## Dec 1978      1036.1790 809.4404 1262.9176 689.41229 1382.9458
## Jan 1979      1036.5608 798.9077 1274.2140 673.10172 1400.0200
## Feb 1979      1070.7985 822.4674 1319.1295 691.00884 1450.5881
## Mar 1979      1128.1704 869.3688 1386.9721 732.36741 1523.9735
## Apr 1979      1058.8200 789.7308 1327.9092 647.28362 1470.3564
## May 1979      1027.8601 748.6463 1307.0739 600.83943 1454.8808
## Jun 1979      1096.1687 806.9760 1385.3613 653.88668 1538.4507
## Jul 1979       948.4537 649.4133 1247.4940 491.11096 1405.7963
## Aug 1979       540.5015 231.7322  849.2708  68.27953 1012.7235
## Sep 1979       984.2214 665.8308 1302.6121 497.28500 1471.1579
## Oct 1979      1073.2842 745.3706 1401.1978 571.78355 1574.7848
## Nov 1979      1010.9422 673.5955 1348.2890 495.01498 1526.8695
## Dec 1979      1078.6251 731.9279 1425.3222 548.39750 1608.8526

ffclog

##          Point Forecast    Lo 80     Hi 80    Lo 95     Hi 95
## Jan 1978      1001.5787 917.1031 1093.8356 875.3079 1146.0653
## Feb 1978      1043.3341 920.1458 1183.0148 860.9358 1264.3755
## Mar 1978      1117.0709 956.5437 1304.5378 881.1262 1416.1959
## Apr 1978      1028.6811 858.7304 1232.2666 780.4440 1355.8755
## May 1978       991.8158 809.2026 1215.6395 726.5660 1353.9014
## Jun 1978      1075.0175 858.7239 1345.7907 762.4430 1515.7365
## Jul 1978       887.9740 695.3665 1133.9313 610.9448 1290.6202
## Aug 1978       376.0420 288.9758  489.3406 251.3710  562.5456
## Sep 1978       935.5515 706.0783 1239.6028 608.3538 1438.7297
## Oct 1978      1047.5181 776.9480 1412.3134 663.2788 1654.3483
## Nov 1978       971.4874 708.5182 1332.0587 599.4900 1574.3178
## Dec 1978      1052.7117 755.2809 1467.2712 633.5387 1749.2253
## Jan 1979      1061.0231 749.1736 1502.6824 623.1212 1806.6629
## Feb 1979      1105.2567 768.2982 1589.9975 633.7602 1927.5307
## Mar 1979      1183.3697 810.0803 1728.6730 662.8238 2112.7242
## Apr 1979      1089.7340 734.8276 1616.0527 596.4747 1990.8977
## May 1979      1050.6807 698.0663 1581.4113 562.2049 1963.5722
## Jun 1979      1138.8205 745.6497 1739.3047 595.8999 2176.3927
## Jul 1979       940.6758 607.0950 1457.5494 481.4816 1837.8084
## Aug 1979       398.3604 253.4576  626.1046 199.5046  795.4251
## Sep 1979       991.0771 621.7528 1579.7819 485.7649 2022.0352
## Oct 1979      1109.6889 686.5230 1793.6900 532.4238 2312.8371
## Nov 1979      1029.1458 627.9594 1686.6395 483.4573 2190.7647
## Dec 1979      1115.1908 671.2080 1852.8541 513.0216 2424.1678

resnaive<-as.numeric(residuals(ffcna))
resdrift<-as.numeric(residuals(ffcdrf))
reslog<-as.numeric(residuals(ffclog, type="response"))

resnaive

##   [1]          NA   14.087976    4.260481    6.058267   -3.986977
##   [6]   41.276687  -21.297725   71.795227  -43.758775  -39.201226
##  [11]  -14.822136   14.266251   42.689355  -17.508147    7.264877
##  [16]  -12.915638   13.854441    7.891326  -26.834853   76.750824
##  [21]  -46.571198  -29.877854   52.247474  -34.397367   21.508452
##  [26]  -10.251114  -11.538570   53.462672   51.115074  -63.852261
##  [31]   69.786793  -31.735426  -13.872467   44.104068   35.295635
##  [36]  -99.467842    4.893686   50.231413   -3.821463   -1.329382
##  [41]  -12.509964   42.607826  -48.026649    2.278019   -9.860960
##  [46]   -2.686969   45.046978  -16.760156   24.903987    4.090877
##  [51]  -26.222418   34.269208  -21.239359  -25.262281   78.672714
##  [56]  -66.882389   35.796694  -52.632497   25.840831   51.241804
##  [61]    9.005797  -31.567817   39.396002  -14.380400  -23.529331
##  [66]   18.403836   18.011408  -13.611350   20.091039   58.544718
##  [71]   -7.311049  -52.115940   34.605660  -23.689264   37.255671
##  [76]  -62.610108   65.114596  -26.903435  -66.797286  116.210163
##  [81]  -58.294141   20.560254  -58.087608  114.898456 -122.828567
##  [86]   49.046452   26.606612  -51.639827    1.654514  -15.974538
##  [91]  -21.307390   58.402863   36.258172   -4.907849  -19.529139
##  [96]  -28.618859   86.913588  -69.292252  -17.435865   55.013441
## [101]   -6.533582   36.272414   48.463562 -123.182776   53.795147
## [106]   23.329583  -54.722134   37.004346   -3.084105   87.243552
## [111]  -69.998825   64.780591  -62.524904   70.405604  -23.562815
## [116]  -74.960675   63.758220  -21.182576   20.754110   60.841352

resdrift

##   [1]           NA   10.5508069    0.7233116    2.5210978   -7.5241465
##   [6]   37.7395174  -24.8348942   68.2580572  -47.2959444  -42.7383959
##  [11]  -18.3593055   10.7290816   39.1521853  -21.0453165    3.7277077
##  [16]  -16.4528079   10.3172711    4.3541561  -30.3720228   73.2136542
##  [21]  -50.1083670  -33.4150237   48.7103049  -37.9345369   17.9712821
##  [26]  -13.7882834  -15.0757396   49.9255022   47.5779044  -67.3894306
##  [31]   66.2496233  -35.2725959  -17.4096366   40.5668984   31.7584652
##  [36] -103.0050114    1.3565161   46.6942432   -7.3586328   -4.8665514
##  [41]  -16.0471338   39.0706569  -51.5638185   -1.2591504  -13.3981291
##  [46]   -6.2241385   41.5098089  -20.2973259   21.3668174    0.5537080
##  [51]  -29.7595879   30.7320386  -24.7765285  -28.7994508   75.1355447
##  [56]  -70.4195589   32.2595245  -56.1696661   22.3036618   47.7046345
##  [61]    5.4686272  -35.1049867   35.8588327  -17.9175694  -27.0665008
##  [66]   14.8666666   14.4742389  -17.1485198   16.5538700   55.0075488
##  [71]  -10.8482182  -55.6531098   31.0684901  -27.2264331   33.7185015
##  [76]  -66.1472776   61.5774267  -30.4406047  -70.3344555  112.6729934
##  [81]  -61.8313104   17.0230844  -61.6247775  111.3612864 -126.3657363
##  [86]   45.5092822   23.0694430  -55.1769965   -1.8826558  -19.5117078
##  [91]  -24.8445594   54.8656935   32.7210030   -8.4450189  -23.0663085
##  [96]  -32.1560280   83.3764187  -72.8294215  -20.9730342   51.4762716
## [101]  -10.0707515   32.7352447   44.9263923 -126.7199451   50.2579778
## [106]   19.7924140  -58.2593033   33.4671765   -6.6212744   83.7063822
## [111]  -73.5359941   61.2434216  -66.0620738   66.8684348  -27.0999842
## [116]  -78.4978448   60.2210502  -24.7197451   17.2169407   57.3041827

reslog

##   [1]           NA   15.1708134   11.9429160   -8.0326480   -9.0840077
##   [6]   45.4983598  -37.6389197    8.2550547    5.0362736  -30.9521327
##  [11]  -23.5287900   15.7354144   41.2554565  -18.9710195   11.6583563
##  [16]  -23.6238210    8.9447242    9.4865022  -42.6149127   11.8956336
##  [21]   -1.7476194  -21.4830768   43.2016529  -33.7186645   18.5882808
##  [26]  -12.0057968   -8.4045137   42.0934920   47.0699967  -69.3422323
##  [31]   59.0813818  -24.6817928    8.2394918   47.1994907   29.9016144
##  [36] -100.8329936    0.5769391   47.9562273   -7.6879115   -5.0935811
##  [41]  -16.2144622   40.1310473  -45.4908876  -17.9301898   32.0386610
##  [46]    1.3593707   36.9419306  -16.5496483   18.6547351   -0.1396053
##  [51]  -30.8815644   29.1242055  -24.6680009  -28.2922168   70.7120498
##  [56]  -44.2118457   70.5625702  -54.9122887   17.0954396   52.4345775
##  [61]    1.8262950  -36.2903358   34.9616929  -14.0344622  -27.4241238
##  [66]   16.1959216   18.5150448    0.3197144   -3.9951770   52.7580566
##  [71]   -9.3977468  -55.2895180   26.8650061  -29.1087647   29.9082167
##  [76]  -59.6017596   60.4266677  -33.7079782  -62.3821652   94.7685793
##  [81] -186.3642523   15.1982770  -61.9757834  116.9782966 -132.8690685
##  [86]   46.6057537   19.0752689  -49.4186475   -3.1497116  -17.8415758
##  [91]  -26.0529387   33.0813322    5.2298801  -14.0713945  -22.7430318
##  [96]  -29.7496747   78.0266077  -76.1067746  -24.4296540   53.1053157
## [101]   -9.5298066   30.7299638   60.7653454  -51.7731881   43.7467877
## [106]   13.6092560  -55.2736092   35.7903037  -12.8760129   81.4646322
## [111]  -86.3411683   69.1839575  -62.3064004   62.0537183    0.3683485
## [116]   -6.7666694  -10.2546999  -37.3282934   21.2053091   52.3527669

resnaive<-na.omit(resnaive)
resdrift<-na.omit(resdrift)
reslog<-na.omit(reslog)

RMSEnaive<-sqrt(sum((resnaive)^2)/119)
RMSEdrift<-sqrt(sum((resdrift)^2)/119)
RMSElog<-sqrt(sum((reslog)^2)/119)

RMSEnaive

## [1] 46.6923

RMSEdrift

## [1] 46.55813

RMSElog

## [1] 45.38429

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.