library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.8
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   2.1.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(fable)
## Loading required package: fabletools
library(tsibble)
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(feasts)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:tsibble':
## 
##     interval
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(tsibbledata)
library(fpp3)
library(ggplot2)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
NYNQGSP_data <- read.csv("//Users//alisonpan//Desktop//Predictive analysis//NYNQGSP.csv", header=TRUE)

##auto.arima model VS ETS model

NYNQGSP_ts <- ts(NYNQGSP_data$NYNQGSP,start=c(2005,1),end=c(2021,7), frequency=12)
plot(NYNQGSP_ts)

train_data = ts(NYNQGSP_data$NYNQGSP,start=c(2005,1), end=c(2021,7), frequency=12)
test_data = ts(NYNQGSP_data$NYNQGSP,start=c(2021,8), frequency=12)

##ARIMA MODEL

ARIMA_model = auto.arima(train_data)
ARIMA_forecast = forecast(ARIMA_model, h = 6)
summary(ARIMA_model)
## Series: train_data 
## ARIMA(0,1,0) 
## 
## sigma^2 = 9.159e+09:  log likelihood = -2551.82
## AIC=5105.63   AICc=5105.65   BIC=5108.92
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE   MAPE      MASE       ACF1
## Training set 4130.082 95464.26 31794.24 0.01931578 2.5052 0.1417862 -0.0247163

Accuracy of ARIMA MODEL

ARIMA_accuracy = accuracy(ARIMA_forecast, test_data)
print(ARIMA_accuracy)
##                       ME      RMSE       MAE          MPE     MAPE      MASE
## Training set    4130.082  95464.26  31794.24   0.01931578  2.50520 0.1417862
## Test set     -780452.133 780969.06 780452.13 -75.66325459 75.66325 3.4804208
##                    ACF1 Theil's U
## Training set -0.0247163        NA
## Test set      0.5472033  45.92973

##ETS MODEL

ETS_model = ets(train_data)
ETS_forecast = forecast(ETS_model, h = 6)
summary(ETS_model)
## ETS(M,N,N) 
## 
## Call:
##  ets(y = train_data) 
## 
##   Smoothing parameters:
##     alpha = 0.9865 
## 
##   Initial states:
##     l = 1061590.7529 
## 
##   sigma:  0.0524
## 
##      AIC     AICc      BIC 
## 5504.942 5505.065 5514.822 
## 
## Training set error measures:
##                    ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 3829.034 95565.75 32132.75 -0.01606958 2.534682 0.1432958
##                     ACF1
## Training set -0.01251638

Accuracy of ETS MODEL

ETS_accuracy = accuracy(ETS_forecast, test_data)
print(ETS_accuracy)
##                       ME      RMSE       MAE          MPE      MAPE      MASE
## Training set    3829.034  95565.75  32132.75  -0.01606958  2.534682 0.1432958
## Test set     -779953.096 780470.35 779953.10 -75.61492237 75.614922 3.4781953
##                     ACF1 Theil's U
## Training set -0.01251638        NA
## Test set      0.54720326  45.90013

##hold out 6 months of data,the auto.arima forecasting VS ETS forecasting

#auto.arima forecasting

ARIMA_forecast
##          Point Forecast   Lo 80   Hi 80   Lo 95   Hi 95
## Aug 2021        1813748 1691097 1936399 1626169 2001326
## Sep 2021        1813748 1640293 1987202 1548472 2079024
## Oct 2021        1813748 1601310 2026185 1488852 2138643
## Nov 2021        1813748 1568446 2059050 1438591 2188905
## Dec 2021        1813748 1539492 2088004 1394310 2233186
## Jan 2022        1813748 1513316 2114180 1354276 2273219
autoplot(NYNQGSP_ts) +
  autolayer(ARIMA_forecast, series = "ARIMA FORECAST")

##ETS forecasting

autoplot(NYNQGSP_ts) +
  autolayer(ETS_forecast, series = "ETS FORECAST")