Exercise 7.8
Zeng Yao
4/27/2021
# read the data
pigs## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1980 76378 71947 33873 96428 105084 95741 110647 100331 94133 103055
## 1981 76889 81291 91643 96228 102736 100264 103491 97027 95240 91680
## 1982 76892 85773 95210 93771 98202 97906 100306 94089 102680 77919
## 1983 81225 88357 106175 91922 104114 109959 97880 105386 96479 97580
## 1984 90974 98981 107188 94177 115097 113696 114532 120110 93607 110925
## 1985 103069 103351 111331 106161 111590 99447 101987 85333 86970 100561
## 1986 82719 79498 74846 73819 77029 78446 86978 75878 69571 75722
## 1987 63292 59380 78332 72381 55971 69750 85472 70133 79125 85805
## 1988 69069 79556 88174 66698 72258 73445 76131 86082 75443 73969
## 1989 66269 73776 80034 70694 81823 75640 75540 82229 75345 77034
## 1990 75982 78074 77588 84100 97966 89051 93503 84747 74531 91900
## 1991 81022 78265 77271 85043 95418 79568 103283 95770 91297 101244
## 1992 93866 95171 100183 103926 102643 108387 97077 90901 90336 88732
## 1993 73292 78943 94399 92937 90130 91055 106062 103560 104075 101783
## 1994 82413 83534 109011 96499 102430 103002 91815 99067 110067 101599
## 1995 88905 89936 106723 84307 114896 106749 87892 100506
## Nov Dec
## 1980 90595 101457
## 1981 101259 109564
## 1982 93561 117062
## 1983 109490 110191
## 1984 103312 120184
## 1985 89543 89265
## 1986 64182 77357
## 1987 81778 86852
## 1988 78139 78646
## 1989 78589 79769
## 1990 81635 89797
## 1991 114525 101139
## 1992 83759 99267
## 1993 93791 102313
## 1994 97646 104930
## 1995#forecasting using the ses() function
pig_model <- ses(pigs , h= 4)
summary(pig_model)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = pigs, h = 4)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249 0.01282239
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.8#95% prediction interval
standard_deviation <- sd(pig_model$residuals)
standard_deviation## [1] 10273.69print(paste0("Lower Confidence Interval: ", pig_model$mean[1]-1.96*standard_deviation))## [1] "Lower Confidence Interval: 78679.9672534162"print(paste0("Upper Confidence Interval: ", pig_model$mean[1]+1.96*standard_deviation))## [1] "Upper Confidence Interval: 118952.844969765"#The forecast of the next observation in the series
autoplot(pig_model) +
autolayer(fitted(pig_model), series = "Fitted")+
ylab("Number of Pigs Slaughtered")
# Question 3
# function takes argument as a vector, because optim() requires vector as the first argument of the function.
sse_fn <- function( pars=c(alpha,l0), y){
# assigning the initial sse value
sse <- 0
# assigning first parameter to alpha
alpha <- pars[1]
# assigning second parametr to l0
l0 <- pars[2]
y_hat <- l0
# calculating sse
for(index in 1:length(y)){
sse <- sse + (y[index] - y_hat)^2
y_hat <- alpha*y[index] + (1 - alpha)*y_hat
}
return(sse)
}
# finding the optimal values of a and l0 using the optim(). Setting initial value of alpha = 0, l0 = first observation value of pigs data set
opt_par<- optim(par = c(0, pigs[1]), y = pigs, fn = sse_fn)
opt_par## $par
## [1] 0.297086 77272.075070
##
## $value
## [1] 19765613447
##
## $counts
## function gradient
## 179 NA
##
## $convergence
## [1] 0
##
## $message
## NULL#Question 4
# finding the optimal values of a and l0 using the optim()
pig_model$model$par## alpha l
## 2.971488e-01 7.726006e+04#Question 5 ##a
#loading the data books
books## Time Series:
## Start = 1
## End = 30
## Frequency = 1
## Paperback Hardcover
## 1 199 139
## 2 172 128
## 3 111 172
## 4 209 139
## 5 161 191
## 6 119 168
## 7 195 170
## 8 195 145
## 9 131 184
## 10 183 135
## 11 143 218
## 12 141 198
## 13 168 230
## 14 201 222
## 15 155 206
## 16 243 240
## 17 225 189
## 18 167 222
## 19 237 158
## 20 202 178
## 21 186 217
## 22 176 261
## 23 232 238
## 24 195 240
## 25 190 214
## 26 182 200
## 27 222 201
## 28 217 283
## 29 188 220
## 30 247 259head(books)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## Paperback Hardcover
## 1 199 139
## 2 172 128
## 3 111 172
## 4 209 139
## 5 161 191
## 6 119 168autoplot(books) +
ggtitle('Sales of Books at a Store') +
xlab('Day') +
ylab('Books Sold')
##b
sesfitp <- ses(books[,1])
sesfith <- ses(books[,2])
summary(sesfitp)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 1])
##
## Smoothing parameters:
## alpha = 0.1685
##
## Initial states:
## l = 170.8271
##
## sigma: 34.8183
##
## AIC AICc BIC
## 318.9747 319.8978 323.1783
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303 -0.2117522
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 207.1097 162.4882 251.7311 138.8670 275.3523
## 32 207.1097 161.8589 252.3604 137.9046 276.3147
## 33 207.1097 161.2382 252.9811 136.9554 277.2639
## 34 207.1097 160.6259 253.5935 136.0188 278.2005
## 35 207.1097 160.0215 254.1979 135.0945 279.1249
## 36 207.1097 159.4247 254.7946 134.1818 280.0375
## 37 207.1097 158.8353 255.3840 133.2804 280.9389
## 38 207.1097 158.2531 255.9663 132.3899 281.8294
## 39 207.1097 157.6777 256.5417 131.5099 282.7094
## 40 207.1097 157.1089 257.1105 130.6400 283.5793autoplot(sesfitp) +
autolayer(fitted(sesfitp), series='Fitted') +
ggtitle('SES Fit and Forecast of Paperback Sales') +
xlab('Day') +
ylab('Books Sale')
summary(sesfith)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, 2])
##
## Smoothing parameters:
## alpha = 0.3283
##
## Initial states:
## l = 149.2861
##
## sigma: 33.0517
##
## AIC AICc BIC
## 315.8506 316.7737 320.0542
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887 -0.1417763
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 239.5601 197.2026 281.9176 174.7799 304.3403
## 32 239.5601 194.9788 284.1414 171.3788 307.7414
## 33 239.5601 192.8607 286.2595 168.1396 310.9806
## 34 239.5601 190.8347 288.2855 165.0410 314.0792
## 35 239.5601 188.8895 290.2306 162.0662 317.0540
## 36 239.5601 187.0164 292.1038 159.2014 319.9188
## 37 239.5601 185.2077 293.9124 156.4353 322.6848
## 38 239.5601 183.4574 295.6628 153.7584 325.3618
## 39 239.5601 181.7600 297.3602 151.1625 327.9577
## 40 239.5601 180.1111 299.0091 148.6406 330.4795autoplot(sesfith) +
autolayer(fitted(sesfith), series='Fitted') +
ggtitle('SES Fit and Forecast of Hardcover Sales') +
xlab('Day') +
ylab('Books Sale')
##c
#Question 6 ##a
#a
holtfitp <- holt(books[,1], h=4)
forecast(holtfitp)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 209.4668 166.6035 252.3301 143.9130 275.0205
## 32 210.7177 167.8544 253.5811 145.1640 276.2715
## 33 211.9687 169.1054 254.8320 146.4149 277.5225
## 34 213.2197 170.3564 256.0830 147.6659 278.7735holtfith <- holt(books[,2], h=4)
forecast(holtfith)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 31 250.1739 212.7390 287.6087 192.9222 307.4256
## 32 253.4765 216.0416 290.9113 196.2248 310.7282
## 33 256.7791 219.3442 294.2140 199.5274 314.0308
## 34 260.0817 222.6468 297.5166 202.8300 317.3334##b
#b##c
#c
sesfitp <- ses(books[,1], h=4)
sesfith <- ses(books[,1], h=4)
autoplot(books[,1]) +
autolayer(holtfitp, series='Holts Method', PI=F) +
autolayer(sesfitp, series='Simple ETS', PI=F) +
ggtitle('Paperback Sales') +
xlab('Day') +
ylab('Books Sales') +
guides(colour=guide_legend(title="Forecast"))
autoplot(books[,2]) +
autolayer(holtfith, series='Holts Method', PI=F) +
autolayer(sesfith, series='Simple ETS', PI=F) +
ggtitle('Hardcover Sales') +
xlab('Day') +
ylab('Books Sales') +
guides(colour=guide_legend(title="Forecast"))
##d
#d
rmsep <- 31.14
ptholtp <- 209.4668
ptsesp <- 207.1097
lowerp <- ptholtp - 1.96 * rmsep
upperp <- ptholtp + 1.96 * rmsep
holtlowerp <- 143.9130
holtupperp <- 275.0205
seslowerp <- 138.8670
sesupperp <- 275.3523
rmseh <- 27.19
ptholth <- 250.1739
ptsesh <- 239.5601
lowerh <- ptholth - 1.96 * rmseh
upperh <- ptholth + 1.96 * rmseh
holtlowerh <- 192.9222
holtupperh <- 307.4256
seslowerh <- 174.7799
sesupperh <- 304.3403
df <- data.frame(c(ptholtp, lowerp, upperp), c(ptholtp, holtlowerp, holtupperp), c(ptsesp, seslowerp, sesupperp), c(ptholth, lowerh, upperh), c(ptholth, holtlowerh, holtupperh), c(ptsesh, seslowerh, sesupperh))
df[4,] <- df[3,] - df[2,]
colnames(df) <- c('Calculated', 'R - holt', 'R - ses', 'Calculated', 'R - holt', 'R - ses')
row.names(df) <- c('Point Forecast', 'Lower 95%', 'Upper 95%', 'Interval Range')
kable(df) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
add_header_above(c(' ', 'Paperback Forecast' = 3, 'Hardcover Forecost' = 3))| Calculated | R - holt | R - ses | Calculated | R - holt | R - ses | |
|---|---|---|---|---|---|---|
| Point Forecast | 209.4668 | 209.4668 | 207.1097 | 250.1739 | 250.1739 | 239.5601 |
| Lower 95% | 148.4324 | 143.9130 | 138.8670 | 196.8815 | 192.9222 | 174.7799 |
| Upper 95% | 270.5012 | 275.0205 | 275.3523 | 303.4663 | 307.4256 | 304.3403 |
| Interval Range | 122.0688 | 131.1075 | 136.4853 | 106.5848 | 114.5034 | 129.5604 |
#Question 7
default <- holt(eggs, h=100)
damped <- holt(eggs, h=100, damped = T)
exponential <- holt(eggs, h=100, exponential = T)
lambda <- holt(eggs, h=100, lambda = 'auto', biasadj = T)
da_ex <- holt(eggs, h=100, exponential = T, damped = T)
da_la <- holt(eggs, h=100, damped = T, lambda = 'auto', biasadj = T)
autoplot(eggs) +
autolayer(default, series='Default', PI=F) +
autolayer(damped, series='Damped', PI=F) +
autolayer(exponential, series='Exponential', PI=F) +
autolayer(lambda, series='Box-Cox Transformed', PI=F) +
autolayer(da_ex, series='Damped & Exponential', PI=F) +
autolayer(da_la, series='Damped & Box-Cox', PI=F) +
ggtitle('Forecast of US Eggs Prices') +
xlab('Year') +
ylab('Price of Dozen Eggs')