About Data Analysis Report

This RMarkdown file contains the report of the data analysis done for the project on forecasting daily bike rental demand using time series models in R. It contains analysis such as data exploration, summary statistics and building the time series models. The final report was completed on Sun Oct 22 02:22:30 2023.

Data Description:

This dataset contains the daily count of rental bike transactions between years 2011 and 2012 in Capital bikeshare system with the corresponding weather and seasonal information.

Data Source: https://archive.ics.uci.edu/ml/datasets/bike+sharing+dataset

Relevant Paper:

Fanaee-T, Hadi, and Gama, Joao, ‘Event labeling combining ensemble detectors and background knowledge’, Progress in Artificial Intelligence (2013): pp. 1-15, Springer Berlin Heidelberg

Task One: Load and explore the data

Load data and install packages

## Import required packages
library("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("ggplot2")
library("dbplyr")
## 
## Attaching package: 'dbplyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     ident, sql
library("timetk")
library("tseries")
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library("forecast")

Describe and explore the data

day_data <- read.csv("D:\\day.csv")
day_data[,"dteday"] <- as.Date(day_data[,"dteday"])
day_data[,"ncnt"] <- day_data[,"cnt"] / max(day_data[,"cnt"])
day_data[,"nr"] <- day_data[,"registered"] / max(day_data[,"registered"])
day_data[,"rr"] <- day_data[,"cnt"] / max(day_data[,"registered"])
summary(day_data)
##     instant          dteday               season            yr        
##  Min.   :  1.0   Min.   :2011-01-01   Min.   :1.000   Min.   :0.0000  
##  1st Qu.:183.5   1st Qu.:2011-07-02   1st Qu.:2.000   1st Qu.:0.0000  
##  Median :366.0   Median :2012-01-01   Median :3.000   Median :1.0000  
##  Mean   :366.0   Mean   :2012-01-01   Mean   :2.497   Mean   :0.5007  
##  3rd Qu.:548.5   3rd Qu.:2012-07-01   3rd Qu.:3.000   3rd Qu.:1.0000  
##  Max.   :731.0   Max.   :2012-12-31   Max.   :4.000   Max.   :1.0000  
##       mnth          holiday           weekday        workingday   
##  Min.   : 1.00   Min.   :0.00000   Min.   :0.000   Min.   :0.000  
##  1st Qu.: 4.00   1st Qu.:0.00000   1st Qu.:1.000   1st Qu.:0.000  
##  Median : 7.00   Median :0.00000   Median :3.000   Median :1.000  
##  Mean   : 6.52   Mean   :0.02873   Mean   :2.997   Mean   :0.684  
##  3rd Qu.:10.00   3rd Qu.:0.00000   3rd Qu.:5.000   3rd Qu.:1.000  
##  Max.   :12.00   Max.   :1.00000   Max.   :6.000   Max.   :1.000  
##    weathersit         temp             atemp              hum        
##  Min.   :1.000   Min.   :0.05913   Min.   :0.07907   Min.   :0.0000  
##  1st Qu.:1.000   1st Qu.:0.33708   1st Qu.:0.33784   1st Qu.:0.5200  
##  Median :1.000   Median :0.49833   Median :0.48673   Median :0.6267  
##  Mean   :1.395   Mean   :0.49538   Mean   :0.47435   Mean   :0.6279  
##  3rd Qu.:2.000   3rd Qu.:0.65542   3rd Qu.:0.60860   3rd Qu.:0.7302  
##  Max.   :3.000   Max.   :0.86167   Max.   :0.84090   Max.   :0.9725  
##    windspeed           casual         registered        cnt      
##  Min.   :0.02239   Min.   :   2.0   Min.   :  20   Min.   :  22  
##  1st Qu.:0.13495   1st Qu.: 315.5   1st Qu.:2497   1st Qu.:3152  
##  Median :0.18097   Median : 713.0   Median :3662   Median :4548  
##  Mean   :0.19049   Mean   : 848.2   Mean   :3656   Mean   :4504  
##  3rd Qu.:0.23321   3rd Qu.:1096.0   3rd Qu.:4776   3rd Qu.:5956  
##  Max.   :0.50746   Max.   :3410.0   Max.   :6946   Max.   :8714  
##       ncnt                nr                 rr          
##  Min.   :0.002525   Min.   :0.002879   Min.   :0.003167  
##  1st Qu.:0.361717   1st Qu.:0.359488   1st Qu.:0.453786  
##  Median :0.521919   Median :0.527210   Median :0.654765  
##  Mean   :0.516909   Mean   :0.526371   Mean   :0.648481  
##  3rd Qu.:0.683498   3rd Qu.:0.687662   3rd Qu.:0.857472  
##  Max.   :1.000000   Max.   :1.000000   Max.   :1.254535

Task Two: Create interactive time series plots

## Read about the timetk package
# ?timetk
#day_data
day_data %>% group_by(yr) %>% 
  plot_time_series(dteday, temp, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Temperature", 
                              .title = "Normalized Temperature vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% 
  plot_time_series(dteday, hum, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Humidity", 
                   .title = "Normalized Humidity vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% 
  plot_time_series(dteday, windspeed, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Windspeed", 
                              .title = "Windspeed vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% 
  plot_time_series(dteday, ncnt, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Bike Rentals", 
                              .title = "Normalized Bike Rentals vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% plot_time_series(dteday, nr, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Registered Rentals", 
                              .title = "Normalized Registered Rentals vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% plot_time_series(dteday, rr, .color_var = season, .x_lab = "Date", .y_lab = "Ratio", 
                              .title = "Ratio of Rentals to Registration vs Date for Day Data", .interactive = TRUE)

Task Three: Smooth time series data

since the humidity and wind speed data has so much noise, I will focus on temperature and bike rentals

data cleaning

day_data[,"temp"] <- tsclean(day_data[,"temp"])
day_data[,"ncnt"] <- tsclean(day_data[,"ncnt"])
day_data[,"nr"] <- tsclean(day_data[,"nr"])
day_data[,"rr"] <- tsclean(day_data[,"rr"])
head(day_data)
##   instant     dteday season yr mnth holiday weekday workingday weathersit
## 1       1 2011-01-01      1  0    1       0       6          0          2
## 2       2 2011-01-02      1  0    1       0       0          0          2
## 3       3 2011-01-03      1  0    1       0       1          1          1
## 4       4 2011-01-04      1  0    1       0       2          1          1
## 5       5 2011-01-05      1  0    1       0       3          1          1
## 6       6 2011-01-06      1  0    1       0       4          1          1
##       temp    atemp      hum windspeed casual registered  cnt       ncnt
## 1 0.344167 0.363625 0.805833 0.1604460    331        654  985 0.11303649
## 2 0.363478 0.353739 0.696087 0.2485390    131        670  801 0.09192105
## 3 0.196364 0.189405 0.437273 0.2483090    120       1229 1349 0.15480835
## 4 0.200000 0.212122 0.590435 0.1602960    108       1454 1562 0.17925178
## 5 0.226957 0.229270 0.436957 0.1869000     82       1518 1600 0.18361258
## 6 0.204348 0.233209 0.518261 0.0895652     88       1518 1606 0.18430112
##           nr        rr
## 1 0.09415491 0.1418082
## 2 0.09645839 0.1153182
## 3 0.17693637 0.1942125
## 4 0.20932911 0.2248776
## 5 0.21854305 0.2303484
## 6 0.21854305 0.2312122
#day_data
day_data %>% group_by(yr) %>% plot_time_series(dteday, temp, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Temperature", .title = "Normalized Temperature vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% plot_time_series(dteday, ncnt, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Bike Rentals", .title = "Normalized Bike Rentals vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% plot_time_series(dteday, nr, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Registered Rentals", .title = "Normalized Registered Rentals vs Date for Day Data", .interactive = TRUE)
day_data %>% group_by(yr) %>% plot_time_series(dteday, rr, .color_var = season, .x_lab = "Date", .y_lab = "Normalized Registered Rentals", .title = "Ratio of Rentals to Registration vs Date for Day Data", .interactive = TRUE)

Task Four: Decompse and access the stationarity of time series data

#day data
day_data[,"temp"] %>% adf.test()
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -1.6785, Lag order = 9, p-value = 0.7144
## alternative hypothesis: stationary
day_data[,"ncnt"] %>% adf.test()
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -1.3084, Lag order = 9, p-value = 0.871
## alternative hypothesis: stationary
day_data[,"nr"] %>% adf.test()
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -2.418, Lag order = 9, p-value = 0.4014
## alternative hypothesis: stationary
day_data[,"rr"] %>% adf.test()
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -1.3084, Lag order = 9, p-value = 0.871
## alternative hypothesis: stationary
#decomposes the data
freq <- 365

#normalized day rentals
norm_rentals <- ts(day_data[, "nr"], frequency = freq)
decomped1 <- stl(norm_rentals, "periodic")
plot(decomped1$time.series[,2], ylab = "Stationary of the Normalized Rental Reservations", 
     xlab = "Day of the Year")

checkresiduals(decomped1$time.series[, 3])

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 1997.3, df = 146, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 146
#normalized day counts
norm_cnt <- ts(day_data[, "ncnt"], frequency = freq)
decomped2 <- stl(norm_cnt, "periodic")
plot(decomped2$time.series[,2], ylab = "Stationary of the Normalized Rental Counts", 
     xlab = "Day of the Year")

checkresiduals(decomped2$time.series[, 3])

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 1437.8, df = 146, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 146
#normalized day rental rates
norm_rr <- ts(day_data[, "rr"], frequency = freq)
decomped3 <- stl(norm_rr, "periodic")
plot(decomped3$time.series[,2], ylab = "Stationary of the Normalized Rental Counts to Reservations", 
     xlab = "Day of the Year")

checkresiduals(decomped3$time.series[, 3])

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 1437.8, df = 146, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 146
print("-------Noramlized Rental Reservations")
## [1] "-------Noramlized Rental Reservations"
shapiro.test(decomped1$time.series[, 3])
## 
##  Shapiro-Wilk normality test
## 
## data:  decomped1$time.series[, 3]
## W = 0.99554, p-value = 0.03334
print("-------Normalized Count of Rentals")
## [1] "-------Normalized Count of Rentals"
shapiro.test(decomped2$time.series[, 3])
## 
##  Shapiro-Wilk normality test
## 
## data:  decomped2$time.series[, 3]
## W = 0.99702, p-value = 0.1993
print("-------Normalized Ratio of Rentals to Reservations")
## [1] "-------Normalized Ratio of Rentals to Reservations"
shapiro.test(decomped3$time.series[, 3])
## 
##  Shapiro-Wilk normality test
## 
## data:  decomped3$time.series[, 3]
## W = 0.99702, p-value = 0.1993

Task Five: Fit and forecast time series data using ARIMA models

#bike count predictions
fit1 <- auto.arima(norm_cnt, seasonal = TRUE, )
hist(fit1$residuals, xlab = "Residual", ylab = "Distribution", main = "Histogram of Model Errors - Bike Count")

shapiro.test(fit1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.83122, p-value < 2.2e-16
prediction1 <- forecast(fit1, 25)
plot(prediction1, xlab = "Date", ylab = "Normalized Count of Rentals", main = "Prediction of Bike Rental Counts")

#predictions of number of reservations
fit2 <- auto.arima(norm_rentals, seasonal = TRUE, )
hist(fit2$residuals, xlab = "Residual", ylab = "Distribution", main = "Histogram of Model Errors - Rental Count")

shapiro.test(fit2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit2$residuals
## W = 0.85305, p-value < 2.2e-16
prediction2 <- forecast(fit2, 25)
plot(prediction2, xlab = "Date", ylab = "Normalized Rentals", main = "Prediction of Bike Rentals")

#bike count predictions
fit2 <- auto.arima(norm_cnt, seasonal = TRUE, )
hist(fit2$residuals, xlab = "Residual", ylab = "Distribution", main = "Histogram of Model Errors - Count to Rental Ratio")

shapiro.test(fit2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit2$residuals
## W = 0.83122, p-value < 2.2e-16
prediction3 <- forecast(fit2, 25)
plot(prediction3, xlab = "Date", ylab = "Normalized Rental Ratio", main = "Prediction of Bike Rentals to Reservations")

Task Six: Findings and Conclusions

After processing the raw data and using the ARIMA package to model ride share data, I was able to make predictions for the 25 days beyond the current data set. Qualitatively the data shows that was the weather gets warmer the number of bike rentals increase, and over the course of two years the number of rentals increases over the number of rentals from the previous year. As the data terminates at the end of one cycle, I expect the number of rentals to increase to a level higher than it was a year before, which is what the models are predicting. Therefore the results were what I expected - the data appears to oscillate up and down over a 1 year period with the overall data moving towards higher rental numbers.