MeDaScIn2016Hyndman.pdf

Slide 1

Slide 1 text

Rob J Hyndman Automatic FoRecasting using R

Slide 2

Slide 2 text

Outline 1 Forecasting the PBS 2 Exponential smoothing 3 ARIMA models 4 TBATS models 5 Hierarchical time series Automatic FoRecasting using R Forecasting the PBS 2

Slide 3

Slide 3 text

Forecasting the PBS Automatic FoRecasting using R Forecasting the PBS 3

Slide 4

Slide 4 text

Forecasting the PBS The Pharmaceutical Beneﬁts Scheme (PBS) is the Australian government drugs subsidy scheme. Many drugs bought from pharmacies are subsidised to allow more equitable access to modern drugs. The cost to government is determined by the number and types of drugs purchased. Currently nearly 1% of GDP ($14 billion). The total cost is budgeted based on forecasts of drug usage. Automatic FoRecasting using R Forecasting the PBS 4

Slide 5

Slide 5 text

Slide 6

Slide 6 text

Slide 7

Slide 7 text

Slide 8

Slide 8 text

Forecasting the PBS In 2001: $4.5 billion budget, under-forecasted by $800 million. Thousands of products. Seasonal demand. Subject to covert marketing, volatile products, uncontrollable expenditure. Although monthly data available for 10 years, data are aggregated to annual values, and only the ﬁrst three years are used in estimating the forecasts. All forecasts being done with the FORECAST function in MS-Excel! Automatic FoRecasting using R Forecasting the PBS 5

Slide 9

Slide 9 text

Slide 10

Slide 10 text

Slide 11

Slide 11 text

Slide 12

Slide 12 text

Slide 13

Slide 13 text

PBS data Automatic FoRecasting using R Forecasting the PBS 6 Total cost: A03 concession safety net group Time $ thousands 1995 2000 2005 0 200 400 600 800 1000 1200

Slide 14

Slide 14 text

PBS data Automatic FoRecasting using R Forecasting the PBS 6 Total cost: A05 general copayments group Time $ thousands 1995 2000 2005 0 50 100 150 200

Slide 15

Slide 15 text

PBS data Automatic FoRecasting using R Forecasting the PBS 6 Total cost: D01 general copayments group Time $ thousands 1995 2000 2005 0 100 200 300 400 500 600 700

Slide 16

Slide 16 text

PBS data Automatic FoRecasting using R Forecasting the PBS 6 Total cost: S01 general copayments group Time $ thousands 1995 2000 2005 0 500 1000 1500

Slide 17

Slide 17 text

PBS data Automatic FoRecasting using R Forecasting the PBS 6 Total cost: R03 general copayments group Time $ thousands 1995 2000 2005 1000 2000 3000 4000 5000

Slide 18

Slide 18 text

Outline 1 Forecasting the PBS 2 Exponential smoothing 3 ARIMA models 4 TBATS models 5 Hierarchical time series Automatic FoRecasting using R Exponential smoothing 7

Slide 19

Slide 19 text

Slide 20

Slide 20 text

Slide 21

Slide 21 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M N,N: Simple exponential smoothing A,N: Holt’s linear method Automatic FoRecasting using R Exponential smoothing 8

Slide 22

Slide 22 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M N,N: Simple exponential smoothing A,N: Holt’s linear method Ad,N: Additive damped trend method Automatic FoRecasting using R Exponential smoothing 8

Slide 23

Slide 23 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M N,N: Simple exponential smoothing A,N: Holt’s linear method Ad,N: Additive damped trend method A,A: Additive Holt-Winters’ method Automatic FoRecasting using R Exponential smoothing 8

Slide 24

Slide 24 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M N,N: Simple exponential smoothing A,N: Holt’s linear method Ad,N: Additive damped trend method A,A: Additive Holt-Winters’ method A,M: Multiplicative Holt-Winters’ method Automatic FoRecasting using R Exponential smoothing 8

Slide 25

Slide 25 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M There are 9 separate exponential smoothing methods. Automatic FoRecasting using R Exponential smoothing 8

Slide 26

Slide 26 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M There are 9 separate exponential smoothing methods. Each can have an additive or multiplicative error, giving 18 separate models. Automatic FoRecasting using R Exponential smoothing 8

Slide 27

Slide 27 text

Exponential smoothing methods Seasonal Component Trend N A M Component (None) (Additive) (Multiplicative) N (None) N,N N,A N,M A (Additive) A,N A,A A,M Ad (Additive damped) Ad,N Ad,A Ad,M There are 9 separate exponential smoothing methods. Each can have an additive or multiplicative error, giving 18 separate models. Only 15 models are numerically stable. Automatic FoRecasting using R Exponential smoothing 8

Slide 28

Slide 28 text

Cross-validation Traditional evaluation Automatic FoRecasting using R Exponential smoothing 9 q q q q q q q q q q q q q q q q q q q q q q q q q time Training data Test data

Slide 29

Slide 29 text

Cross-validation Traditional evaluation Standard cross-validation Automatic FoRecasting using R Exponential smoothing 9 q q q q q q q q q q q q q q q q q q q q q q q q q time Training data Test data q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Slide 30

Slide 30 text

Cross-validation Traditional evaluation Standard cross-validation Time series cross-validation Automatic FoRecasting using R Exponential smoothing 9 q q q q q q q q q q q q q q q q q q q q q q q q q time Training data Test data q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Slide 31

Slide 31 text

Slide 32

Slide 32 text

Slide 33

Slide 33 text

Akaike’s Information Criterion AICC = −2 log(L) + 2k + 2(k + 1)(k + 2) T − k where L is the likelihood, k is the number of estimated parameters in the model and T is the number of observations in the series. Minimizing the Gaussian AIC is asymptotically equivalent (as T → ∞) to minimizing MSE from one-step forecasts on test set via time series cross-validation. Automatic FoRecasting using R Exponential smoothing 10

Slide 34

Slide 34 text

Slide 35

Slide 35 text

ets algorithm in R Automatic FoRecasting using R Exponential smoothing 11 Based on Hyndman, Koehler, Snyder & Grose (IJF 2002): Apply each of 15 models that are appropriate to the data. Optimize parameters and initial values using MLE. Select best method using AICc. Produce forecasts using best method. Obtain prediction intervals using underlying state space model.

Slide 36

Slide 36 text

Slide 37

Slide 37 text

Slide 38

Slide 38 text

Slide 39

Slide 39 text

Exponential smoothing Automatic FoRecasting using R Exponential smoothing 12 300 400 500 600 1960 1980 2000 Year Number of sheep (millions) level 80 95 Forecasts from ETS(M,A,N)

Slide 40

Slide 40 text

Exponential smoothing Automatic FoRecasting using R Exponential smoothing 13 300 400 500 600 1960 1980 2000 Year Number of sheep (millions) level 80 95 Forecasts from ETS(M,A,N) fit <- ets(livestock) fcast <- forecast(fit) plot(fcast)

Slide 41

Slide 41 text

Exponential smoothing Automatic FoRecasting using R Exponential smoothing 14 0.3 0.6 0.9 1.2 1.5 1995 2000 2005 2010 Year millions of scripts level 80 95 Forecasts from ETS(M,Ad,M)

Slide 42

Slide 42 text

Exponential smoothing Automatic FoRecasting using R Exponential smoothing 15 0.3 0.6 0.9 1.2 1.5 1995 2000 2005 2010 Year millions of scripts level 80 95 Forecasts from ETS(M,Ad,M) fit <- ets(h02) fcast <- forecast(fit) plot(fcast)

Slide 43

Slide 43 text

Outline 1 Forecasting the PBS 2 Exponential smoothing 3 ARIMA models 4 TBATS models 5 Hierarchical time series Automatic FoRecasting using R ARIMA models 16

Slide 44

Slide 44 text

auto.arima algorithm in R Based on Hyndman and Khandakar (JSS 2008): Select no. differences via unit root tests. Use stepwise search to traverse model space, starting with a simple model and considering nearby variants. For each model, optimize parameters using MLE. Select best method using AICc. Produce forecasts and prediction intervals using best method. Automatic FoRecasting using R ARIMA models 17

Slide 45

Slide 45 text

Slide 46

Slide 46 text

Slide 47

Slide 47 text

Slide 48

Slide 48 text

Slide 49

Slide 49 text

Auto ARIMA Automatic FoRecasting using R ARIMA models 18 300 400 500 1960 1980 2000 Year Number of sheep (millions) level 80 95 Forecasts from ARIMA(0,1,0) with drift

Slide 50

Slide 50 text

Auto ARIMA Automatic FoRecasting using R ARIMA models 19 300 400 500 1960 1980 2000 Year Number of sheep (millions) level 80 95 Forecasts from ARIMA(0,1,0) with drift fit <- auto.arima(livestock) fcast <- forecast(fit) plot(fcast)

Slide 51

Slide 51 text

Auto ARIMA Automatic FoRecasting using R ARIMA models 20 0.3 0.6 0.9 1.2 1995 2000 2005 2010 Year millions of scripts level 80 95 Forecasts from ARIMA(3,1,3)(0,1,1)[12]

Slide 52

Slide 52 text

Auto ARIMA Automatic FoRecasting using R ARIMA models 21 0.3 0.6 0.9 1.2 1995 2000 2005 2010 Year millions of scripts level 80 95 Forecasts from ARIMA(3,1,3)(0,1,1)[12] fit <- auto.arima(h02) fcast <- forecast(fit) plot(fcast)

Slide 53

Slide 53 text

Outline 1 Forecasting the PBS 2 Exponential smoothing 3 ARIMA models 4 TBATS models 5 Hierarchical time series Automatic FoRecasting using R TBATS models 22

Slide 54

Slide 54 text

TBATS model TBATS Trigonometric terms for seasonality Box-Cox transformations for heterogeneity ARMA errors for short-term dynamics Trend (possibly damped) Seasonal (including multiple and non-integer periods) Automatic algorithm described in De Livera, Hyndman and Snyder (JASA 2011). Automatic FoRecasting using R TBATS models 23

Slide 55

Slide 55 text

Examples Automatic FoRecasting using R TBATS models 24 7 8 9 10 1995 2000 2005 Year Thousands of barrels level 80 95 Forecasts from TBATS(1, {0,0}, 1, {<52.1785714285714,9>}) fit <- tbats(gas) fcast <- forecast(fit) plot(fcast)

Slide 56

Slide 56 text

Examples Automatic FoRecasting using R TBATS models 25 0 200 400 0 10 20 30 Weeks Number of calls level 80 95 Forecasts from TBATS(0.607, {0,0}, −, {<169,5>, <845,4>}) fit <- tbats(callcentre) fcast <- forecast(fit) plot(fcast)

Slide 57

Slide 57 text

Outline 1 Forecasting the PBS 2 Exponential smoothing 3 ARIMA models 4 TBATS models 5 Hierarchical time series Automatic FoRecasting using R Hierarchical time series 26

Slide 58

Slide 58 text

Hierarchical time series A hierarchical time series is a collection of several time series that are linked together in a hierarchical structure. Total A AA AB AC B BA BB BC C CA CB CC Examples Manufacturing product hierarchies Sales by state and region Automatic FoRecasting using R Hierarchical time series 27

Slide 59

Slide 59 text

Slide 60

Slide 60 text

Slide 61

Slide 61 text

Optimal reconciliation Traditional approaches: Bottom-up Top-down Middle-out A better way: 1 Forecast all series at all levels of aggregation. 2 Reconcile the forecasts by making the smallest possible changes such that they add up. 3 Extremely fast algorithm implemented in the hts package for R. Automatic FoRecasting using R Hierarchical time series 28

Slide 62

Slide 62 text

Slide 63

Slide 63 text

Slide 64

Slide 64 text

Slide 65

Slide 65 text

Slide 66

Slide 66 text

Slide 67

Slide 67 text

Slide 68

Slide 68 text

Australian tourism Automatic FoRecasting using R Hierarchical time series 29

Slide 69

Slide 69 text

Australian tourism Automatic FoRecasting using R Hierarchical time series 29 Hierarchy: States (7) Zones (27) Regions (82)

Slide 70

Slide 70 text

Australian tourism Automatic FoRecasting using R Hierarchical time series 29 Hierarchy: States (7) Zones (27) Regions (82) Base forecasts ETS (exponential smoothing) models

Slide 71

Slide 71 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: Total Year Visitor nights 1998 2000 2002 2004 2006 2008 60000 65000 70000 75000 80000 85000

Slide 72

Slide 72 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: NSW Year Visitor nights 1998 2000 2002 2004 2006 2008 18000 22000 26000 30000

Slide 73

Slide 73 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: VIC Year Visitor nights 1998 2000 2002 2004 2006 2008 10000 12000 14000 16000 18000

Slide 74

Slide 74 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: Nth.Coast.NSW Year Visitor nights 1998 2000 2002 2004 2006 2008 5000 6000 7000 8000 9000

Slide 75

Slide 75 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: Metro.QLD Year Visitor nights 1998 2000 2002 2004 2006 2008 8000 9000 11000 13000

Slide 76

Slide 76 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: Sth.WA Year Visitor nights 1998 2000 2002 2004 2006 2008 400 600 800 1000 1200 1400

Slide 77

Slide 77 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: X201.Melbourne Year Visitor nights 1998 2000 2002 2004 2006 2008 4000 4500 5000 5500 6000

Slide 78

Slide 78 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: X402.Murraylands Year Visitor nights 1998 2000 2002 2004 2006 2008 0 100 200 300

Slide 79

Slide 79 text

Base forecasts Automatic FoRecasting using R Hierarchical time series 30 Domestic tourism forecasts: X809.Daly Year Visitor nights 1998 2000 2002 2004 2006 2008 0 20 40 60 80 100

Slide 80

Slide 80 text

Reconciled forecasts Automatic FoRecasting using R Hierarchical time series 31 Total 2000 2005 2010 65000 80000 95000

Slide 81

Slide 81 text

Reconciled forecasts Automatic FoRecasting using R Hierarchical time series 31 NSW 2000 2005 2010 18000 24000 30000 VIC 2000 2005 2010 10000 14000 18000 QLD 2000 2005 2010 14000 20000 Other 2000 2005 2010 18000 24000

Slide 82

Slide 82 text

Reconciled forecasts Automatic FoRecasting using R Hierarchical time series 31 Sydney 2000 2005 2010 4000 7000 Other NSW 2000 2005 2010 14000 22000 Melbourne 2000 2005 2010 4000 5000 Other VIC 2000 2005 2010 6000 12000 GC and Brisbane 2000 2005 2010 6000 9000 Other QLD 2000 2005 2010 6000 12000 Capital cities 2000 2005 2010 14000 20000 Other 2000 2005 2010 5500 7500

Slide 83

Slide 83 text

Forecast evaluation Select models using all observations; Re-estimate models using ﬁrst 12 observations and generate 1- to 8-step-ahead forecasts; Increase sample size one observation at a time, re-estimate models, generate forecasts until the end of the sample; In total 24 1-step-ahead, 23 2-steps-ahead, up to 17 8-steps-ahead for forecast evaluation. Automatic FoRecasting using R Hierarchical time series 32

Slide 84

Slide 84 text

Slide 85

Slide 85 text

Slide 86

Slide 86 text

Slide 87

Slide 87 text

Hierarchy: states, zones, regions MAPE h = 1 h = 2 h = 4 h = 6 h = 8 Average Top Level: Australia Bottom-up 3.79 3.58 4.01 4.55 4.24 4.06 OLS 3.83 3.66 3.88 4.19 4.25 3.94 WLS 3.68 3.56 3.97 4.57 4.25 4.04 Level: States Bottom-up 10.70 10.52 10.85 11.46 11.27 11.03 OLS 11.07 10.58 11.13 11.62 12.21 11.35 WLS 10.44 10.17 10.47 10.97 10.98 10.67 Level: Zones Bottom-up 14.99 14.97 14.98 15.69 15.65 15.32 OLS 15.16 15.06 15.27 15.74 16.15 15.48 WLS 14.63 14.62 14.68 15.17 15.25 14.94 Bottom Level: Regions Bottom-up 33.12 32.54 32.26 33.74 33.96 33.18 OLS 35.89 33.86 34.26 36.06 37.49 35.43 WLS 31.68 31.22 31.08 32.41 32.77 31.89 Automatic FoRecasting using R Hierarchical time series 33

Slide 88

Slide 88 text

Forecast package history Pre 2003 Collection of functions used for consulting projects July/August 2003 ets and thetaf added August 2006 v1.0 available on CRAN May 2007 auto.arima added July 2008 JSS paper (Hyndman & Khandakar) September 2009 v2.0. Unbundled. May 2010 arﬁma added Feb/March 2011 tslm, stlf, naive, snaive added August 2011 v3.0. Box Cox transformations added December 2011 tbats added April 2012 Package moved to github November 2012 v4.0. nnetar added June 2013 Major speed-up of ets January 2014 v5.0. tsoutliers and tsclean added May 2015 v6.0. Added several new plots December 2015 264,000 package downloads in one month! February 2016 v7.0. Added ggplot2 graphics & bias adjustment Automatic FoRecasting using R Hierarchical time series 34

Slide 89

Slide 89 text

For further information robjhyndman.com Slides for this talk. Links to all papers and books. Links to R packages. A blog about forecasting research. OTexts.org/fpp Free online book based on forecast package for R. Automatic FoRecasting using R Hierarchical time series 35