library(readxl)
library(tidyquant)
library(fpp3)
library(moments)
library(tsibble)
library(tsibbledata)
library(ggplot2)
library(dplyr)
library(seasonal)
library(seasonalview)
library(fabletools)
library(fable)
library(ggfortify)
library(quantmod)
library(GGally)
library(feasts)
raw <- read_excel("C:/Users/harle/OneDrive/Desktop/Class Files/Fall 2022/Forecasting/TSLA MAX Years.xlsx", 
                             col_types = c("date", "numeric"))

teslamod <- raw %>% 
  mutate(row=row_number()) %>% 
  as_tsibble(index='row') %>%
  slice(2648:3120) #sliced the data to only included 2021 and 2022 since those values matter much more according to Nathan, the finance guy.

ggplot(teslamod, mapping = aes(Date, Close)) +
  geom_line() +
  labs(title = "Tesla Stock Closing Price", subtitle = "Initial Plot")

teslamod %>%
  gg_tsdisplay(difference(Close), plot_type='partial')
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

teslarima <- teslamod %>%
  model(arima010 = ARIMA(Close ~ pdq(0,2,0)),
        stepwise = ARIMA(Close),
        search = ARIMA(Close, stepwise=FALSE))
#Differentiated twice to remove any trend. The partial acf only showed significant lags far into the past so I set p = 0. I did not do a moving average so q = 0.
glance(teslarima) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 3 × 6
##   .model   sigma2 log_lik   AIC  AICc   BIC
##   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 search     108.  -1771. 3556. 3556. 3585.
## 2 stepwise   110.  -1779. 3559. 3559. 3563.
## 3 arima010   230.  -1949. 3900. 3900. 3904.
teslarima %>% select(search) %>% gg_tsresiduals()

#The residuals look really good. They resemble white noise and are normally distributed.
restest <- teslamod %>%
  features(difference(Close), ljung_box, lag = 10)
restest
## # A tibble: 1 × 2
##   lb_stat lb_pvalue
##     <dbl>     <dbl>
## 1    19.9    0.0299
teslacast <- teslarima %>%
  forecast(h=1) %>%
  filter(.model=='search')
teslacast %>%
  hilo(80)
## # A tsibble: 1 x 5 [1]
## # Key:       .model [1]
##   .model   row       Close .mean                  `80%`
##   <chr>  <dbl>      <dist> <dbl>                 <hilo>
## 1 search  3121 N(186, 108)  186. [172.7612, 199.3468]80
teslacast %>%
  autoplot(teslamod)

#This forecast looks good and looks similar to the NAIVE method. 

teslanaive <- teslamod %>%
  model('Simple' = ETS(Close ~ error("A") + trend("N") + season("N")))

teslanaivecast <- teslanaive %>%
  forecast(h = 1)
teslanaivecast %>%
  hilo(80)
## # A tsibble: 1 x 5 [1]
## # Key:       .model [1]
##   .model   row       Close .mean                `80%`
##   <chr>  <dbl>      <dist> <dbl>               <hilo>
## 1 Simple  3121 N(187, 110)  187. [173.82, 200.6799]80
teslanaivecast %>%
  autoplot(teslamod)

#The ARIMA and NAIVE methods gave the exact same point forecast of 187. I'm going with the 80% prediction intervals of the ARIMA of [172.7612, 199.3468].