Exponential Smoothing
Use the ses() function in R to find the optimal values of α and ℓ 0, and generate forecasts for the next four months.
alpha = 0.2971 l = 77260.0561
# Simple exponential smoothing model is used to generate the next four months.
fc <- ses(pigs, h = 4)
sumfc <- summary(fc)##
## 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
# Plot of forecast and fitted values in the original time series below to provide visualization.
autoplot(pigs) + autolayer(fitted(fc), series = "Fitted values") + autolayer(fc, series = "Forecast", PI = FALSE) + labs(x = "Year", y = "")Compute a 95% prediction interval for the first forecast using ^ y ± 1.96 s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
σ = 10308.58
# Computation of residual standard error (rse)
s1 <- sqrt(sum(residuals(fc)^2) / (length(residuals(fc)) - 2))
# rse from model
s2 <- sqrt(fc$model$sigma2)
c("S1" = s1, "S2" = s2)## S1 S2
## 10308.58 10308.58
These figures agrees with the residual standard error as output in the model summary above. Next we use ^ y ± 1.96 s for the 95% prediction interval.
In this case, the first prediction interval is slightly different from the 95% prediction interval in the model summary indicated above because a precises estimate was not used. When a precise estimate is used for the critical value from the normal distribution (1.959964), the prediction interval is the same as to the output in the model summary.
## Lower95 Upper95
## 78611.6 119021.2
## Lower95 Upper95
## 78611.97 119020.84
Plot the series and discuss the main features of the data.
Use the ses() function to forecast each series, and plot the forecasts.
#SES function to forecast each series.
fc1 <- ses(books[ , "Paperback"], h = 4)
fc2 <- ses(books[ , "Hardcover"], h = 4)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, "Paperback"], h = 4)
##
## 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
##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = books[, "Hardcover"], h = 4)
##
## 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
#Plot of forecast
autoplot(books) + autolayer(fc1, PI = FALSE, series = "Paperback") + autolayer(fc2, PI = FALSE, series = "Hardcover") + labs(x = "Day", y = "")Compute the RMSE values for the training data in each case.
#Compute RMSE training data
rmse1 <- sqrt(mean(residuals(fc1)^2))
rmse2 <- sqrt(mean(residuals(fc2)^2))
c(Paperback = rmse1, Hardcover = rmse2)## Paperback Hardcover
## 33.63769 31.93101
To verify these figures, the accuracy function can be used. The RMSE when using the accuracy function matches the function above.
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303 -0.2117522
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887 -0.1417763
Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
#Holt's linear method
fc3 <- holt(books[ , "Paperback"], h = 4)
fc4 <- holt(books[ , "Hardcover"], h = 4)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = books[, "Paperback"], h = 4)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 170.699
## b = 1.2621
##
## sigma: 33.4464
##
## AIC AICc BIC
## 318.3396 320.8396 325.3456
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.717178 31.13692 26.18083 -5.508526 15.58354 0.6602122
## ACF1
## Training set -0.1750792
##
## Forecasts:
## 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.7735
##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = books[, "Hardcover"], h = 4)
##
## Smoothing parameters:
## alpha = 1e-04
## beta = 1e-04
##
## Initial states:
## l = 147.7935
## b = 3.303
##
## sigma: 29.2106
##
## AIC AICc BIC
## 310.2148 312.7148 317.2208
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1357882 27.19358 23.15557 -2.114792 12.1626 0.6908555
## ACF1
## Training set -0.03245186
##
## Forecasts:
## 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
#Comparison plot
autoplot(books) + autolayer(fc1, PI = FALSE, series = "Paperback - Simple") + autolayer(fc2, PI = FALSE, series = "Hardcover - Simple") + autolayer(fc3, PI = FALSE, series = "Paperback - Holt") + autolayer(fc4, PI = FALSE, series = "Hardcover - Holt") + labs(x = "Day", y = "")From the visiualization, the Holt series shows an upwards trend as opposed to the Simple.
Compare the RMSE measures of Holt’s method for the two series to those of simple exponential smoothing in the previous question. (Remember that Holt’s method is using one more parameter than SES.) Discuss the merits of the two forecasting methods for these data sets.
When comparing the RMSE measures of Holt’s method to the two series of simple exponential, it can be determined that the linear model produces lower training RMSEs for both time series. Both time series show a consistent upward trend, and the linear model includes a parameter to account for this trend. It can be expected that the linear model would result in a more accurate forecast when outside the training data.
rmse3 <- sqrt(mean(residuals(fc3)^2))
rmse4 <- sqrt(mean(residuals(fc4)^2))
c(Paperback = rmse3, Hardcover = rmse4)## Paperback Hardcover
## 31.13692 27.19358
## ME RMSE MAE MPE MAPE MASE
## Training set -3.717178 31.13692 26.18083 -5.508526 15.58354 0.6602122
## ACF1
## Training set -0.1750792
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1357882 27.19358 23.15557 -2.114792 12.1626 0.6908555
## ACF1
## Training set -0.03245186
#The first column represent Paperback and the second column represents Hardcover.
rbind(Simple = c(rmse1, rmse2), Holt = c(rmse3, rmse4)) ## [,1] [,2]
## Simple 33.63769 31.93101
## Holt 31.13692 27.19358
Compare the forecasts for the two series using both methods. Which do you think is best?
When comparing the two series, Holt’s linear model is preferable since both time series exhibit a greater upward trend.
#Forecast Comparison
cbind(Day = 1:4,
"Paper-Simple" = fc1$mean,
"Paper-Holt" = fc3$mean,
"Hard-Simple" = fc2$mean,
"Hard-Holt" = fc4$mean)## Time Series:
## Start = 31
## End = 34
## Frequency = 1
## Day Paper-Simple Paper-Holt Hard-Simple Hard-Holt
## 31 1 207.1097 209.4668 239.5601 250.1739
## 32 2 207.1097 210.7177 239.5601 253.4765
## 33 3 207.1097 211.9687 239.5601 256.7791
## 34 4 207.1097 213.2197 239.5601 260.0817
Calculate a 95% prediction interval for the first forecast for each series, using the RMSE values and assuming normal errors. Compare your intervals with those produced using ses and holt.
Below is the comparison of the prediction intervals using the RMSE values which assuming normally distributed errors compared to the prediction intervals as reported in the model summaries using ses and holt.
There are slight differences between model output and the estimates using RMSE. When RSE is used instead of the RMSE the prediction intervals match the values from the model output.
# Method using approximation rmse
approx1 <- fc1$mean[1] + c(-1, 1) * 1.96 * rmse1
approx2 <- fc2$mean[1] + c(-1, 1) * 1.96 * rmse2
approx3 <- fc3$mean[1] + c(-1, 1) * 1.96 * rmse3
approx4 <- fc4$mean[1] + c(-1, 1) * 1.96 * rmse4
# Method using precise RSE
est1 <- fc1$mean[1] + c(-1, 1) * qnorm(0.975) * sqrt(fc1$model$sigma2)
est2 <- fc2$mean[1] + c(-1, 1) * qnorm(0.975) * sqrt(fc2$model$sigma2)
est3 <- fc3$mean[1] + c(-1, 1) * qnorm(0.975) * sqrt(fc3$model$sigma2)
est4 <- fc4$mean[1] + c(-1, 1) * qnorm(0.975) * sqrt(fc4$model$sigma2)
rbind("Simple-FromModel" = c(fc1$lower[1,2], fc1$upper[1,2], fc2$lower[1,2], fc2$upper[1,2]),
"Simple-FromRMSE" = c(approx1, approx2),
"Simple-FromRSE" = c(est1, est2),
"Holt-FromModel" = c(fc3$lower[1,2], fc3$upper[1,2], fc4$lower[1,2], fc4$upper[1,2]),
"Holt-FromRMSE" = c(approx3, approx4),
"Holt-FromRSE" = c(est3, est4)) ## 95% 95% 95% 95%
## Simple-FromModel 138.8670 275.3523 174.7799 304.3403
## Simple-FromRMSE 141.1798 273.0395 176.9753 302.1449
## Simple-FromRSE 138.8670 275.3523 174.7799 304.3403
## Holt-FromModel 143.9130 275.0205 192.9222 307.4256
## Holt-FromRMSE 148.4384 270.4951 196.8745 303.4733
## Holt-FromRSE 143.9130 275.0205 192.9222 307.4256
Which model gives the best RMSE?
From plotting the time series, it can be determined that there is a strong downward trend, which will be accounted for using Holt’s linear trend model. Since a Box-Cox transformation will be experimented with the linear trend models a transformed time series on the same plot is displayed. The optimal parameter for the Box-Cox transformation is λ=0.3956183. The holt function estimates the optimal choice of the parameters α and β (for all models), as well as ϕ (for the models with damped trend).H=100 was set in order to clearly see the differences in the forecasts produced by the four models. From the plot below, the Holt model should not be used since it is observed that the prices go negative. The basic Holt model forecasts prices indicate with a downward trend which then goes negative in 2016. The BoxCox and BC-bias models generates forecasts that trend downward however does not indicate negative prices.The Damped-BC-bias model generates forecasts that trend upward, therefore it is the model that gives the best RMSE. The linear model with Box-Cox transformation and bias-adjustment performs the best which has an RMSE of 26.39.
## [1] 0.3956183
cbind(Original = eggs, "Box-Cox" = BoxCox(eggs, lambda = "auto")) %>% autoplot(facet = TRUE) + labs( x = "Year", y = "")# basic holt
eggsfc1 <- holt(eggs, h = 100)
# holt box-cox
eggsfc2 <- holt(eggs, h = 100, lambda = "auto", biasadj = FALSE)
# holt box-cox with bias-adjustment
eggsfc2a <- holt(eggs, h = 100, lambda = "auto", biasadj = TRUE)
# holt damped trend
eggsfc3 <- holt(eggs, h = 100, damped = TRUE)
# holt box-cox & damped trend
eggsfc4 <- holt(eggs, h = 100, lambda = "auto", biasadj = FALSE, damped = TRUE)
# holt box-cox with bias-adjustment & damped trend
eggsfc4a <- holt(eggs, h = 100, lambda = "auto", biasadj = TRUE, damped = TRUE)
# plot all forecasts
autoplot(eggs) + autolayer(eggsfc1, series = "Holt", PI = FALSE, lty = 1, lwd = 1.5) + autolayer(eggsfc2, series = "BoxCox", PI = FALSE, lty = 2, lwd = 1.5) + autolayer(eggsfc2a, series = "BC-bias", PI = FALSE, lty = 3, lwd = 1.5) + autolayer(eggsfc3, series = "Damped", PI = FALSE, lty = 1, lwd = 1.5) + autolayer(eggsfc4, series = "Damped-BC", PI = FALSE, lty = 2, lwd = 1.5) + autolayer(eggsfc4a, series = "Damped-BC-bias", PI = FALSE, lty = 3, lwd = 1.5) + labs(x = "Year", y = "") + geom_hline(yintercept = 0, lty = 3, lwd = 1) #BC-bias is the best.
c(Holt = accuracy(eggsfc1)[2],
BoxCox = accuracy(eggsfc2)[2],
"BC-bias" = accuracy(eggsfc2a)[2],
Damped = accuracy(eggsfc3)[2],
"Damped-BC" = accuracy(eggsfc4)[2],
"Damped-BC-bias" = accuracy(eggsfc4a)[2])## Holt BoxCox BC-bias Damped Damped-BC
## 26.58219 26.39376 26.38689 26.54019 26.53321
## Damped-BC-bias
## 26.58589
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
colID <- colnames(retaildata)[8]
myts <- ts(retaildata[ , colID], frequency=12, start=c(1982,4))
autoplot(myts) + labs(x = "Year", y = "")Why is multiplicative seasonality necessary for this series?
Base on the plot above, it appears that multiplicative seasonality will be more effective in modeling the time series due to the variability in the series is increasing over time.
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
Both forecasts are displayed and each of the models produce forecasts that are very close to each other do a good job of replicating the historical seasonality of the time series.
retailfc1 <- hw(myts, seasonal = "multiplicative")
retailfc2 <- hw(myts, seasonal = "multiplicative", damped = TRUE)
autoplot(myts) + autolayer(retailfc1, series = "HW", PI = FALSE) + autolayer(retailfc2, series = "Damped", PI = FALSE) + labs(x = "Year", y = "")Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
The training RMSEs of both models are similar however the damped model RMSE is slightly lower (13.15) than the Holt-Winter(HW) model RMSE (13.17).Across all metrics the damped Holt-Winters model is stronger than the basic Holt-Winters model.
## HW Damped
## 13.17456 13.15306
hw <- c(accuracy(retailfc1)[2], sqrt(retailfc1$model$sigma2),
retailfc1$model$aic, retailfc1$model$aicc, retailfc1$model$bic)
damped <- c(accuracy(retailfc2)[2], sqrt(retailfc2$model$sigma2),
retailfc2$model$aic, retailfc2$model$aicc, retailfc2$model$bic)
metrics <- round(cbind(HW = hw, Damped = damped), 3)
rownames(metrics) <- c("RMSE", "RSE", "AIC", "AICc", "BIC")
metrics## HW Damped
## RMSE 13.175 13.153
## RSE 0.086 0.080
## AIC 4326.553 4266.723
## AICc 4328.239 4268.612
## BIC 4393.580 4337.693
Check that the residuals from the best method look like white noise.
Checkresiduals function is used to check the residuals of the damped Holt-Winters model. The residuals don’t resemble white noise and they exhibit strong auto-correlation with 1, 6, and 12-month lags. The Ljung-Box test confirms that the residuals do not appear to be white noise.
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 143.08, df = 8, p-value < 2.2e-16
##
## Model df: 16. Total lags used: 24
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 150.7, df = 7, p-value < 2.2e-16
##
## Model df: 17. Total lags used: 24
Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 8 in Section 3.7?
The close-up plot compared to the seasonal naive model does worse as well as worse compared to the damped Holt-Winters on the training set. On the test set it can seen that the damped Holt-Winters model under-predicts actuals and does worse than the seasonal naive model.
# Training set
myts.train <- window(myts, end = c(2010, 12))
# HW model
hw.fc <- hw(myts.train, seasonal = "multiplicative", damped = TRUE, h = 36)
hw.rmse <- accuracy(hw.fc, myts)[ , 2]
# Seasonal naive model
snaive.fc <- snaive(myts.train, h=36)
snaive.rmse <- accuracy(snaive.fc, myts)[ , 2]
# RMSE
cbind("Damped HW" = hw.rmse, "Naive" = snaive.rmse)## Damped HW Naive
## Training set 13.14634 26.30758
## Test set 23.10869 20.91121
Considering the RMSE metrics above the damped Holt-Winters model performs much better than the seasonal naive model on the training set (RMSE of 13.15 vs. 26.31) however it performs worse on the test set (23.1 vs. 20.9).
autoplot(myts) + autolayer(hw.fc, series = "Damped-HW", PI = FALSE) + autolayer(snaive.fc, series = "Seas-Naive", PI = FALSE) + autolayer(fitted(hw.fc), series = "DampedHW") + autolayer(fitted(snaive.fc), series = "SeasNaive") + labs(x = "Year", y = "") + geom_vline(xintercept = 2011, lty = 3, lwd = 1) + coord_cartesian(xlim = c(2008, 2013), ylim = c(225, 425))How does that compare with your best previous forecasts on the test set?
Comparing accuracy measures for this model both with and without robust fitting to the models in the previous model. It appears that the seasonal naive model again performs the best on the test set. The version without robust fitting performs better on both the training and test sets. Comparing the STL with ETS models with and without robust fitting to the seasonal naive model, it appears the seasonal naive model is more responsive to changes even though it has a tendency to aim too high.
# STL fit on training set
stl.fc <- stlf(myts.train, h = 36, s.window = 13, lambda = "auto", method = "ets", robust = FALSE)
# Summary
stl.fc$model## ETS(A,N,N)
##
## Call:
## ets(y = na.interp(x), model = etsmodel, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.7218
##
## Initial states:
## l = 68.0462
##
## sigma: 11.1992
##
## AIC AICc BIC
## 3686.948 3687.019 3698.479
# Robust fitting
stl2.fc <- stlf(myts.train, h = 36, s.window = 13, lambda = "auto", method = "ets", robust = TRUE)
stl2.fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2011 318.6545 303.2389 334.0649 295.0761 342.2207
## Feb 2011 295.7317 277.4166 314.0391 267.7177 323.7273
## Mar 2011 314.6487 293.8492 335.4388 282.8345 346.4407
## Apr 2011 302.9312 279.9009 325.9493 267.7042 338.1297
## May 2011 307.4998 282.4438 332.5415 269.1738 345.7925
## Jun 2011 298.6160 271.6783 325.5368 257.4110 339.7815
## Jul 2011 301.0873 272.3972 329.7583 257.2013 344.9288
## Aug 2011 318.2419 287.9113 348.5524 271.8464 364.5904
## Sep 2011 331.8900 300.0025 363.7564 283.1128 380.6174
## Oct 2011 344.8448 311.4736 378.1935 293.7981 395.8389
## Nov 2011 349.9492 315.1532 384.7213 296.7226 403.1195
## Dec 2011 404.4045 368.2744 440.5121 349.1385 459.6179
## Jan 2012 318.6545 281.1384 356.1398 261.2648 375.9721
## Feb 2012 295.7317 256.9163 334.5118 236.3524 355.0279
## Mar 2012 314.6487 274.6138 354.6483 253.4045 375.8099
## Apr 2012 302.9312 261.6842 344.1391 239.8316 365.9391
## May 2012 307.4998 265.0909 349.8680 242.6223 372.2817
## Jun 2012 298.6160 255.0621 342.1257 231.9857 365.1425
## Jul 2012 301.0873 256.4293 345.6992 232.7675 369.2987
## Aug 2012 318.2419 272.5229 363.9152 248.2997 388.0768
## Sep 2012 331.8900 285.1311 378.6032 260.3573 403.3151
## Oct 2012 344.8448 297.0684 392.5751 271.7560 417.8254
## Nov 2012 349.9492 301.1692 398.6820 275.3248 424.4625
## Dec 2012 404.4045 354.6894 454.0770 328.3527 480.3566
## Jan 2013 318.6545 267.8808 369.3719 240.9765 396.2000
## Feb 2013 295.7317 243.9768 347.4236 216.5494 374.7655
## Mar 2013 314.6487 261.9834 367.2527 234.0751 395.0779
## Apr 2013 302.9312 249.3318 356.4645 220.9265 384.7802
## May 2013 307.4998 253.0032 361.9290 224.1225 390.7186
## Jun 2013 298.6160 243.2185 353.9420 213.8584 383.2048
## Jul 2013 301.0873 244.8185 357.2828 214.9964 387.0054
## Aug 2013 318.2419 261.1361 375.2764 230.8720 405.4436
## Sep 2013 331.8900 273.9546 389.7551 243.2520 420.3622
## Oct 2013 344.8448 286.0908 403.5291 254.9553 434.5701
## Nov 2013 349.9492 290.3782 409.4497 258.8098 440.9224
## Dec 2013 404.4045 344.0870 464.6593 312.1279 496.5339
# RMSE
stl.rmse <- accuracy(stl.fc, myts)[ , 2]
stl2.rmse <- accuracy(stl2.fc, myts)[ , 2]
# Comparison
cbind("Damped HW" = hw.rmse, "Naive" = snaive.rmse, "STL/ETS" = stl.rmse, "STL/ETS/Robust" = stl2.rmse)## Damped HW Naive STL/ETS STL/ETS/Robust
## Training set 13.14634 26.30758 10.76922 12.03861
## Test set 23.10869 20.91121 23.97219 27.87913
# Plot
autoplot(myts, lwd = 1) + autolayer(snaive.fc, series = "Seas-Naive", PI = FALSE) + autolayer(stl.fc, series = "STL/ETS", PI = FALSE) + autolayer(stl2.fc, series = "STL/ETS/Robust", PI = FALSE) + labs( x = "Year", y = "") + geom_vline(xintercept = 2011, lty = 3, lwd = 1) + coord_cartesian(xlim = c(2008, 2013), ylim = c(225, 425))