Time Series in R

Time Series and R

Agenda

Our Data
Time Series Objects
Stationarity
Autocorrelation
ARIMA
Prophet
Neural Networks
XGBoost
Accuracy

ts

zoo

xts

Convert to `xts`

citi <- readr::read_csv(file.path(dataDir, 'citibike_201307_201707.csv'))

citiTS <- xts(citi %>% dplyr::select(Trips), order.by=citi$Date)

citiTS['2017-07-23/']

           Trips
2017-07-23 47779
2017-07-24 44702
2017-07-25 66620
2017-07-26 71672
2017-07-27 67066
2017-07-28 65089
2017-07-29 53897
2017-07-30 60402
2017-07-31 65523

Stationarity

\(E[X_t] = E[X_{t+h}]\)
\(\gamma_x(h, 0) = \gamma_x(h+t, t)\)

Autocorrelation Function

\[ \gamma_{x}(h) = Cov(X_{t+h}, X_t) \]

\[ \rho_x(h) = \frac{\gamma_{x}(h)}{\gamma_{x}(0)} = Cor(X_{t+h}, X_t) \]

ACF

acf(citiTS)

Partial Autocorrelation Function

\[ \alpha(h) = \phi_{hh} \] where \(\phi_{hh}\) is the last element of \[ \phi_{h} = \Gamma^{-1}_h(\gamma_h) \] and \[ \Gamma_h = [\gamma(i-j)]_{i,j=1}^h \] and \[ \gamma_h = [\gamma(1), \gamma(2), \dots, \gamma(h)]^T \]

PACF

pacf(citiTS)

Dickey Fuller Test

\[ DF_\tau = \frac{\hat{\gamma}}{SE(\hat{\gamma})} \]

ADF Test

adf.test(citiTS)

    Augmented Dickey-Fuller Test

data:  citiTS
Dickey-Fuller = -2.9189, Lag order = 11, p-value = 0.1893
alternative hypothesis: stationary

KPSS Test

\[ Y_t = d_t + r_t + \epsilon_t \] \[ r_t = r_{t-1} + u_t \] \[ u_t \sim iid(0, \sigma_u^2) \] Test that \[ \sigma_u = 0 \]

KPSS Test

kpss.test(citiTS)

    KPSS Test for Level Stationarity

data:  citiTS
KPSS Level = 5.9245, Truncation lag parameter = 8, p-value = 0.01

Difference the Data

\[ y_t^{'} = y_t - y_{t-1} \]

Diff

citiDiff <- diff(citiTS)
dygraph(citiDiff, elementId='CitiDiff') %>% dyRangeSelector() %>% dyUnzoom()

ACF and PACF

acf(citiDiff[-1])
pacf(citiDiff[-1])

ADF on Diff

adf.test(citiDiff[-1])

    Augmented Dickey-Fuller Test

data:  citiDiff[-1]
Dickey-Fuller = -17.493, Lag order = 11, p-value = 0.01
alternative hypothesis: stationary

KPSS on Diff

kpss.test(citiDiff[-1])

    KPSS Test for Level Stationarity

data:  citiDiff[-1]
KPSS Level = 0.015114, Truncation lag parameter = 8, p-value = 0.1

Box-Cox Transformation

\[ y(\lambda)= \begin{cases} \frac{y^\lambda - 1}{\lambda},& \text{if } \lambda \neq 0\\ log(y), & \text{if } \lambda = 0 \end{cases} \]

Box-Cox

citiBC <- BoxCox(citiTS, lambda=0)
dygraph(citiBC, elementId='CitiBC') %>% dyRangeSelector() %>% dyUnzoom()

ACF and PACF

acf(citiBC)
pacf(citiBC)

Starter Models

Autoregressive

\[ \text{AR}(p) \]

\[ y_t = \beta_0 + \beta_1y_{t-1} + \beta_2y_{t-2} + \dots + \beta_py_{t-p} + \epsilon_t \]

\[ y_t = c + \sum_{i=1}^p \phi_ix_{t-i} + \epsilon_t = c + \sum_{i=1}^p \phi_iB^ix_{t} + \epsilon_t \]

\[ (1 - \phi_1B^1 - \phi_2B^2 - \dots - \phi_pB^p)y_t = \epsilon_t \]

AR(p)

ar14 <- ar(citiTS, order.max=14)
ar14$ar

 [1]  0.48285396  0.03748923  0.01301544  0.01364710  0.05740610
 [6]  0.10362113  0.13929120 -0.01445352  0.02920200 -0.04967610
[11] -0.05221794  0.05979242  0.05699536  0.09126062

Moving Average

\[ \text{MA}(q) \]

\[ y_t = \epsilon_t + \theta_1\epsilon_{t-1} + \theta_2\epsilon_{t-2} + \ldots + \theta_q\epsilon_{t-q} \]

\[ y_t = (1 + \theta_1B^1 + \theta_2B^2 + \ldots + \theta_qB^q)\epsilon_t \]

MA(q)

ma14 <- ma(citiTS, 14)
summary(ma14)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   5337   21823   30914   30836   39158   60702      14

Autoregressive Moving Average

\[ (1 - \phi_1B^1 - \ldots - \phi_pB^p)y_t = (1 + \theta_1B^1 + \ldots + \theta_qB^q)\epsilon_t \]

ARMA(p, q)

arma1_1 <- arma(citiTS, order=c(1, 1))
summary(arma1_1)

Call:
arma(x = citiTS, order = c(1, 1))

Model:
ARMA(1,1)

Residuals:
   Min     1Q Median     3Q    Max 
-47976  -4214   1077   5069  20960 

Coefficient(s):
            Estimate  Std. Error  t value Pr(>|t|)    
ar1         0.990907    0.004651  213.067   <2e-16 ***
ma1        -0.749162    0.033873  -22.117   <2e-16 ***
intercept 309.551314  152.815760    2.026   0.0428 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fit:
sigma^2 estimated as 63353316,  Conditional Sum-of-Squares = 93889663352,  AIC = 30876.34

Autoregressive Integrated Moving Average

\[ (1 - \phi_1B^1 - \ldots - \phi_pB^p)(1-B)^dy_t = (1 + \theta_1B^1 + \ldots + \theta_qB^q)\epsilon_t \]

ARIMA(p, d, q)

arima1_1_1 <- arima(citiTS, order=c(1, 1, 1))
summary(arima1_1_1)

Call:
arima(x = citiTS, order = c(1, 1, 1))

Coefficients:
         ar1      ma1
      0.3834  -0.9016
s.e.  0.0278   0.0108

sigma^2 estimated as 56481996:  log likelihood = -15340.13,  aic = 30686.27

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set 169.4511 7512.918 5554.759 -14.66341 31.23256 0.9318706
                   ACF1
Training set 0.03266496

Seasonal ARIMA

\[ (1 - \phi_pB^p)(1 - \Phi_PB^P)(1 - B)^d(1-B)^Dy_t = (1 + \theta_qB^q)(1 + \Theta_QB^Q)\epsilon_t \]

ARIMA(p, d ,q)(P, D, Q)

arima_season <- arima(citiTS, order=c(1, 1, 1), seasonal=c(12, 2, 12))
summary(arima_season)

Call:
arima(x = citiTS, order = c(1, 1, 1), seasonal = c(12, 2, 12))

Coefficients:
          ar1      ma1     sar1     sar2     sar3    sar4    sar5    sar6
      -0.4859  -0.9998  -0.2511  -0.0837  -0.0011  0.1572  0.0111  0.3046
s.e.   0.2157   0.0035   0.2579   0.0321   0.0160  0.0234  0.0259  0.0089
        sar7    sar8    sar9    sar10   sar11    sar12     sma1     sma2
      0.3297  0.5942  0.2559  -0.5208  0.0967  -0.1624  -0.7959  -0.5989
s.e.  0.1686  0.1876  0.2031   0.1055  0.0883   0.1959   0.0716   0.0857
        sma3     sma4    sma5     sma6     sma7     sma8    sma9   sma10
      0.0765  -0.0006  0.2784  -0.1645  -0.0077  -0.2307  0.3967  0.9158
s.e.  0.0286   0.0494  0.0534   0.0080   0.1215   0.0658  0.0129  0.0936
        sma11    sma12
      -0.8379  -0.0310
s.e.   0.1090   0.0778

sigma^2 estimated as 48239970:  log likelihood = -15222.75,  aic = 30499.5

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set 244.1609 6938.476 4892.245 -11.09788 28.59577 0.8207268
                    ACF1
Training set -0.01560002

Automatic Fitting

citi_arima <- auto.arima(citiTS, 
                         max.p=10, max.d=2, max.q=10, 
                         max.P=30, max.D=2, max.Q=30, 
                         seasonal=TRUE)

ARIMA Forecast

arimaForecast <- forecast(citi_arima, h=30)

autoplot(arimaForecast)

Prophet

\[ y(t) = g(t) + s(t) + h(t) + \epsilon_t \]

Trend

\[ g(t) = \frac{C}{1 + exp \left( -k(t-m) \right) } \]

Seasonality

\[ s(t) = \sum_{n=1}^{N} \left( a_n cos \left( \frac{2\pi{}nt}{P} \right) + b_n sin \left( \frac{2\pi{}nt}{P} \right) \right) \]

Holidays

\[ Z(t) = \left[ \mathbf{1}(t \in D_1), \ldots, \mathbf{1}(t \in D_L) \right] \]

Data For Prophet

citiDF <- tibble::tibble(ds=index(citiTS), y=as.numeric(citiTS))
citiDF

# A tibble: 1,484 x 2
           ds     y
       <date> <dbl>
 1 2013-07-01 16650
 2 2013-07-02 22745
 3 2013-07-03 21864
 4 2013-07-04 22326
 5 2013-07-05 21842
 6 2013-07-06 20467
 7 2013-07-07 20477
 8 2013-07-08 21615
 9 2013-07-09 26641
10 2013-07-10 25732
# ... with 1,474 more rows

Prophet with Defaults

prophet1 <- prophet(citiDF)

Initial log joint probability = -48.3415
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

futureDF1 <- make_future_dataframe(prophet1, periods=30)
preds1 <- predict(prophet1, df=futureDF1)

Prophet Results

plot.prophet(prophet1, preds1, elementID='preds1') %>% 
    dyRangeSelector(dateWindow=as.Date(c('2017-01-01', '2017-10-30')))

Model Components

prophet_plot_components(prophet1, fcst=preds1)

Specify Seasonality

prophet2 <- prophet(citiDF, 
                    yearly.seasonality=TRUE, 
                    weekly.seasonality=TRUE
)

Initial log joint probability = -48.3415
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

futureDF2 <- make_future_dataframe(prophet2, periods=30)
preds2 <- predict(prophet2, df=futureDF2)

Seasonal Components

prophet_plot_components(prophet2, fcst=preds2)

Holiday Data

holidays <- readr::read_csv(file.path(dataDir, 'holidays_2013_2017.csv'))

tail(holidays)

# A tibble: 6 x 2
          ds          holiday
      <date>            <chr>
1 2017-07-04 Independence Day
2 2017-09-04        Labor Day
3 2017-10-09     Columbus Day
4 2017-11-10     Veterans Day
5 2017-11-23 Thanksgiving Day
6 2017-12-25    Christmas Day

Prophet with Holidays

prophet3 <- prophet(citiDF, 
                    yearly.seasonality=TRUE, 
                    weekly.seasonality=TRUE,
                    holidays=holidays)

Initial log joint probability = -48.3415
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

futureDF3 <- make_future_dataframe(prophet3, periods=30)
preds3 <- predict(prophet3, df=futureDF3)

Seasonal Components

prophet_plot_components(prophet3, fcst=preds3)

Neural Network

Simple Neural Network

net1 <- nnetar(citiTS)

netPreds1 <- forecast(net1, h=30)

Neural Net Forecast

autoplot(netPreds1, include=180)

Boosted Trees

XGBoost

\[ \hat{y}_i^{t} = \sum_{k=1}^t f_k(x_i) = \hat{y}_i^{(t-1)} + f_t(x_i) \]

where

\[ \hat{f_t}(x_t)=\sum_{m=1}^M\hat{c}_mI(x_t \in R_m) \]

XGBoost

citiXG1 <- xgbar(as.ts(citiTS), maxlag=14, 
                 trend_method="differencing", 
                 seas_method="fourier")

citiXGPreds1 <- forecast(citiXG1, h=30)

XGBoost Predictions

autoplot(citiXGPreds1, include=180)

Accuracy

forecastAccuracy <- function(y_hat, y)
{
    mean(abs(as.vector(y_hat) - as.vector(y)))
}

results <- tibble::tibble(
    Method=c('ARIMA', 'Prophet', 'Neural Net', 'XGBoost'),
    MAE=c(
        forecastAccuracy(arimaForecast$mean, citiTest_ts),
        forecastAccuracy(preds3 %>% 
                             filter(ds >= '2017-07-31') %>% pull(yhat), citiTest_ts),
        forecastAccuracy(netPreds1$mean, citiTest_ts),
        forecastAccuracy(citiXGPreds1$mean, citiTest_ts)
    )
)

Accuracy

ggplot(results, aes(x=Method, y=MAE)) + geom_bar(stat='identity', width=.01)

What Have We Seen?

Different time series types in R
Interactive Plots
Stationarity
ACF/PACF
Dickey Fuller and KPSS Tests
Differencing
Box-Cox
ARIMA
Prophet
Neural Nets
XGBoost

Going Foward

Use the forecast package
Fit CNNs and RNNs
timetk
tibbletime
Multivariate Models
GARCH Models

Thank You

Jared P. Lander

Chief Data Scientist of Lander Analytics
Author of R for Everyone
Adjunct Professor at Columbia University
Organizer of New York Open Statistical Programming (The R) Meetup
Organizer of New York R Conference

Agenda

Convert to xts

Stationarity

Autocorrelation Function

ACF

Partial Autocorrelation Function

PACF

Dickey Fuller Test

ADF Test

KPSS Test

KPSS Test

Difference the Data

Diff

ACF and PACF

ADF on Diff

KPSS on Diff

Box-Cox Transformation

Box-Cox

ACF and PACF

Autoregressive

AR(p)

Moving Average

MA(q)

Autoregressive Moving Average

ARMA(p, q)

Autoregressive Integrated Moving Average

ARIMA(p, d, q)

Seasonal ARIMA

ARIMA(p, d ,q)(P, D, Q)

Automatic Fitting

ARIMA Forecast

Prophet

Trend

Seasonality

Holidays

Data For Prophet

Prophet with Defaults

Prophet Results

Model Components

Specify Seasonality

Seasonal Components

Holiday Data

Prophet with Holidays

Seasonal Components

Simple Neural Network

Neural Net Forecast

XGBoost

XGBoost

XGBoost Predictions

Accuracy

Accuracy

What Have We Seen?

Going Foward

Jared P. Lander

Packages

Convert to `xts`