<jupyter_start><jupyter_text>Rental bike forecasting using SARIMAX<jupyter_code>library(readr)
library(dplyr)
library(ggplot2)
library(gridExtra)
library(xts)
library(fpp2)
library(DMwR)
library(Metrics)
library(tseries)
options(repr.plot.width=8, repr.plot.height=4)<jupyter_output><empty_output><jupyter_text>1. Data import<jupyter_code>path <- '../data/bike/bike-sharing-processed.csv'
data <- read_csv(path)
head(data)
tail(data)
nrow(data)<jupyter_output><empty_output><jupyter_text>2. Convert to ts object<jupyter_code># define date range as per the data
date_range <- seq(as.Date('2011-01-01'), as.Date('2012-12-31'), by='day')

# no need to put the date column as it will generate as ts index
dayofyear <- as.numeric(format(date_range[1], '%j'))

bike_data <- ts(data[, 2:14], start=c(2011, dayofyear), frequency = 365)
head(bike_data)
autoplot(bike_data[, 12]) +
    ggtitle('Bike count') +
    ylab('count')<jupyter_output><empty_output><jupyter_text>3. Data analysis<jupyter_code>autoplot(diff(bike_data[, 12])) +
    ggtitle('Differenced bike count') +
    ylab('count')<jupyter_output><empty_output><jupyter_text>- The data seems to be I(1)  3. Train test split<jupyter_code># train and test split
# last 2 month will be test set rest training set
# we will drop the explanatory  variables for baseline models
# as they will not be need for forecasting

test_size <- as.numeric(61)
train_size <- nrow(bike_data)- test_size

train <- head(bike_data, train_size)
test <- tail(bike_data, test_size)

autoplot(train[, 12], series = 'train') + 
    autolayer(test[, 12], series = 'test') +
    ggtitle('Rental bike count') +
    ylab('count') + 
    xlab('Year') +
    guides(colour=guide_legend(title="Data"))
tail(train[, 12])
head(test[, 12])
head(train)
# extract exogenous variables
exog_col <- colnames(train)[8:11]

# hurricane sandy to the exogenous list
# exog_col <- c(exog_col, colnames(train)[13])

train_exog <- train[, exog_col]
test_exog <- test[, exog_col]

head(train_exog)
# head(test_exog)<jupyter_output><empty_output><jupyter_text>4. Order evaluation using auto.arima()<jupyter_code>model <- auto.arima(train[, 12], d = 1, D = 1, max.p = 2, max.q = 2, max.P = 2,
          max.Q = 2, max.d = 1, max.D = 1, start.p = 0,
          start.q = 0, start.P = 0, start.Q = 0)
summary(model)<jupyter_output>Series: train[, 12] 
ARIMA(1,1,1)(0,1,0)[365] 

Coefficients:
         ar1      ma1
      0.2636  -0.9276
s.e.  0.0649   0.0286

sigma^2 estimated as 1666471:  log likelihood=-2612.12
AIC=5230.24   AICc=5230.32   BIC=5241.39

Training set error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE
Training set 11.17378 866.6923 415.1966 -47.38094 54.25316 0.1688453
                    ACF1
Training set 0.005671186<jupyter_text>5. Model fitting<jupyter_code># use arima order determined in previous step
# fit <- Arima(train[, 12], order = c(1, 1, 1), seasonal = c(0, 1, 0),
#             include.mean=FALSE, xreg=train_exog)
# saveRDS(fit, '../checkpoints/bike-sharing/R-models/sarimax-v1.rds')
fit <- readRDS('../checkpoints/bike-sharing/R-models/sarimax-v1.rds')
train_pred <- fitted(fit)                            
test_pred <- forecast(fit, h=61, xreg=test_exog)$mean
autoplot(bike_data[, 12], series = 'actual data') +
    autolayer(train_pred, series = 'train prediction') +
    autolayer(test_pred, series = 'test prediction') +
    xlab('Year') +
    ylab('count') +
    ggtitle('Forecasting using SARIMAX')<jupyter_output><empty_output><jupyter_text>6. Model evaluation<jupyter_code>train_rmse <- sqrt(mse(train[, 12], train_pred))
train_mae <- mae(train[, 12], train_pred)
train_nrmse <- train_rmse/sd(train[, 12])

test_rmse <- sqrt(mse(test[, 12], test_pred))
test_mae <- mae(test[, 12], test_pred)
test_nrmse <- test_rmse/sd(test[, 12])

print(paste0('Training NRMSE :', round(train_nrmse, 3)))
print(paste0('Training MAE :', round(train_mae, 3)))

print(paste0('Test NRMSE :', round(test_nrmse, 3)))
print(paste0('Test MAE :', round(test_mae, 3)))
autoplot(test[, 12], series = 'test data') +
    autolayer(test_pred, series = 'test prediction') +
    xlab('Year') +
    ylab('count') +
    ggtitle('Forecasting using SARIMAX on test data') +
    theme(legend.position="bottom", legend.box = "horizontal")<jupyter_output><empty_output><jupyter_text>7. Residual and co-efficient analysis<jupyter_code>summary(fit)
checkresiduals(fit)<jupyter_output>Series: train[, 12] 
Regression with ARIMA(1,1,1)(0,1,0)[365] errors 

Coefficients:
         ar1      ma1      temp      atemp         hum   windspeed
      0.3181  -0.9214  3722.671   841.1244  -3298.9396  -2589.1698
s.e.  0.0791   0.0436  2035.774  2048.4309    373.3004    634.2482

sigma^2 estimated as 1234455:  log likelihood=-2564.39
AIC=5142.78   AICc=5143.16   BIC=5168.8

Training set error measures:
                   ME     RMSE      MAE       MPE    MAPE      MASE        ACF1
Training set 9.790382 740.9833 354.5853 -36.22509 42.2925 0.1441969 -0.01206287

	Ljung-Box test

data:  Residuals from Regression with ARIMA(1,1,1)(0,1,0)[365] errors
Q* = 325.03, df = 128, p-value < 2.2e-16

Model df: 6.   Total lags used: 134
{'fraction_non_alphanumeric': 0.145655546935609, 'fraction_numerical': 0.08320403413498836, 'mean_word_length': 3.977799227799228, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '505869', 'license': 'MIT', 'n_tokens_mistral': 2330, 'n_tokens_neox': 1978, 'n_words': 432, 'path': '02.bike_count/11.bike-count-SARIMAX.ipynb', 'repo_name': 'mkumar73/time-series', 'size': '5156'}