setwd("C:/Users/artem/OneDrive - UNCG/ERM-413")
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
library(ggplot2)
library(tidyverse)
## Warning: package 'purrr' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ── 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(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

Data Exploration

data <- read.csv("MiniProject3Data.csv")
stream_data <- ts(data, start=c(2019, 1), end=c(2023, 12), frequency = 12)
plot(stream_data)

plot(log(stream_data))

decomp_stream <- decompose(log(stream_data), type = "additive")
plot(decomp_stream)

The stream data has an upward trend. Seasonality occurs every year.

Holt-Winters

hw_stream <- HoltWinters(log(stream_data))
hw_stream
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = log(stream_data))
## 
## Smoothing parameters:
##  alpha: 0.6210134
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a    9.98407559
## b    0.03083690
## s1   0.06113446
## s2   0.09033630
## s3   0.12240674
## s4   0.12726561
## s5   0.09249781
## s6   0.04507197
## s7  -0.04934751
## s8  -0.07630076
## s9  -0.18004175
## s10 -0.17483247
## s11 -0.14239864
## s12 -0.07520257
plot(hw_stream)

hw_forecast <- forecast(hw_stream, h = 12)
plot(hw_forecast)

Alpha (level): 0.62 - High, reacts to change Beta (trend): 0 - Low, not constantly gaining or losing subscribers Gamma (seasonality): 1 - High, strong seasonality, regular subscribers

Interactive Plots

dat <- data.frame(
  Time = time(stream_data)[(length(stream_data) - length(hw_stream$fitted[,1]) + 1):length(stream_data)],
  Observed = as.numeric(stream_data)[(length(stream_data) - length(hw_stream$fitted[,1]) + 1):length(stream_data)],
  Fitted = as.numeric(hw_stream$fitted[,1])
)

p <- plot_ly() %>%
  add_trace(x = time(stream_data), y = as.numeric(stream_data), 
            type = 'scatter', mode = 'lines', 
            name = 'Observed', line = list(color = 'blue')) %>%
  add_trace(x = time(hw_stream$fitted), y = as.numeric(hw_stream$fitted[,1]), 
            type = 'scatter', mode = 'lines', 
            name = 'Fitted', line = list(color = 'red'))

forecast_df <- data.frame(
  Time = seq(max(time(hw_stream$fitted)) + 1/frequency(stream_data), 
             by = 1/frequency(stream_data), 
             length.out = 12),
  Forecast = as.numeric(hw_forecast$mean),
  Lower80 = as.numeric(hw_forecast$lower[,1]),
  Upper80 = as.numeric(hw_forecast$upper[,1]),
  Lower95 = as.numeric(hw_forecast$lower[,2]),
  Upper95 = as.numeric(hw_forecast$upper[,2])
)

p <- p %>%
  add_trace(data = forecast_df, x = ~Time, y = ~Forecast, 
            type = 'scatter', mode = 'lines', 
            name = 'Forecast', line = list(color = 'green')) %>%
  add_ribbons(data = forecast_df, x = ~Time, 
              ymin = ~Lower95, ymax = ~Upper95,
              name = '95% Confidence', 
              fillcolor = 'rgba(200, 200, 200, 0.3)',
              line = list(color = 'transparent')) %>%
  add_ribbons(data = forecast_df, x = ~Time, 
              ymin = ~Lower80, ymax = ~Upper80,
              name = '80% Confidence', 
              fillcolor = 'rgba(150, 150, 150, 0.3)',
              line = list(color = 'transparent')) %>%
  layout(title = 'StreamFlix Forecast',
         xaxis = list(title = 'Year'),
         yaxis = list(title = 'Subscribers'),
         hovermode = 'closest')

p
hw_optimal <- HoltWinters(stream_data, optim.start = c(alpha = 0.3, beta = 0.1, gamma = 0.1))
hw_custom <- HoltWinters(stream_data, 
                         alpha = 0.7, 
                         beta = 0.4,
                         gamma = 0.5)
fitted_optimal <- window(stream_data, start=time(hw_optimal$fitted)[1])
rmse_optimal <- sqrt(mean((fitted_optimal - hw_optimal$fitted[,1])^2))

fitted_custom <- window(stream_data, start=time(hw_custom$fitted)[1])
rmse_custom <- sqrt(mean((fitted_custom - hw_custom$fitted[,1])^2))

cat("RMSE for optimal model:", rmse_optimal, "\n")
## RMSE for optimal model: 703.6056
cat("RMSE for custom model:", rmse_custom, "\n")
## RMSE for custom model: 888.4074
best_model <- if(rmse_optimal < rmse_custom) hw_optimal else hw_custom

hw_forecast <- forecast(best_model, h = 12)

lambda <- BoxCox.lambda(stream_data)
hw_bc <- HoltWinters(BoxCox(stream_data, lambda))
hw_bc_forecast <- forecast(hw_bc, h = 12, lambda=lambda)
## Warning in InvBoxCox(pmean, lambda, biasadj, list(level = level, upper = upper,
## : biasadj information not found, defaulting to FALSE.
q <- plot_ly() %>%
  add_trace(x = time(stream_data), y = as.numeric(stream_data), 
            type = 'scatter', mode = 'lines', 
            name = 'Observed', line = list(color = 'blue', width = 1)) %>%
  add_trace(x = time(best_model$fitted), y = as.numeric(best_model$fitted[,1]), 
            type = 'scatter', mode = 'lines', 
            name = 'Fitted', line = list(color = 'red', width = 2)) %>%
  add_trace(x = time(hw_bc_forecast$mean), y = as.numeric(hw_bc_forecast$mean), 
            type = 'scatter', mode = 'lines', 
            name = 'Forecast', line = list(color = 'green', width = 2)) %>%
  add_ribbons(x = time(hw_bc_forecast$mean), 
              ymin = as.numeric(hw_bc_forecast$lower[,1]), 
              ymax = as.numeric(hw_bc_forecast$upper[,1]),
              name = '95% Confidence', 
              fillcolor = 'rgba(150, 150, 150, 0.3)',
              line = list(color = 'transparent')) %>%
  layout(title = 'StreamFlix Forecast',
         xaxis = list(title = 'Year'),
         yaxis = list(title = 'Subscribers'),
         hovermode = 'closest')

decomp <- decompose(stream_data)
q_decomp <- subplot(
  plot_ly(x = time(decomp$x), y = as.numeric(decomp$x), type = 'scatter', mode = 'lines', name = 'Original'),
  plot_ly(x = time(decomp$trend), y = as.numeric(decomp$trend), type = 'scatter', mode = 'lines', name = 'Trend'),
  plot_ly(x = time(decomp$seasonal), y = as.numeric(decomp$seasonal), type = 'scatter', mode = 'lines', name = 'Seasonal'),
  plot_ly(x = time(decomp$random), y = as.numeric(decomp$random), type = 'scatter', mode = 'lines', name = 'Random'),
  nrows = 4, shareX = TRUE, titleY = TRUE
) %>% layout(title = "StreamFlix HW Decomposition")

q
q_decomp

ARMA

adf.test(stream_data) # p-value = 0.01 - Significant but seasonality
## Warning in adf.test(stream_data): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  stream_data
## Dickey-Fuller = -5.9269, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
stream_diff <- diff(stream_data, lag = 12)
adf.test(stream_diff) # p-value = 0.14 - Insignificant but stationary
## 
##  Augmented Dickey-Fuller Test
## 
## data:  stream_diff
## Dickey-Fuller = -3.0876, Lag order = 3, p-value = 0.1398
## alternative hypothesis: stationary
plot(stream_diff)

plot(decompose(stream_diff))

par(mfrow = c(1,2))
acf(stream_diff)
pacf(stream_diff)