Discussion-4.R

library(fpp2)

## Loading required package: ggplot2

## Loading required package: forecast

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

## Loading required package: fma

## Loading required package: expsmooth

library(readr)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

library(readxl)
library(forecast)
library(Metrics)

## Warning: package 'Metrics' was built under R version 3.6.1

## 
## Attaching package: 'Metrics'

## The following object is masked from 'package:forecast':
## 
##     accuracy

GNP <- read_excel("C:/Users/Joseph/Downloads/GNP.xls")

GNP.ts <- ts(GNP$GNP, start = c(1947,1), end = c(2019,1), frequency = 4)
#1) Plot data to identify unusual trends:
plot(GNP.ts)

#2) Need to transform? 
#No unusual trends spotted, significant trend, no box-cox is needed
#3-5) Using auto.arima to determine parameters (p,d,q)
fit1 <- auto.arima(GNP.ts)
#Determines (1,2,2), no seasonal effect
print(paste("AICc:",round(fit1$aicc,2), "BIC:", round(fit1$bic,2)))

## [1] "AICc: 3116.78 BIC: 3131.28"

#ets model
fit2 <- ets(GNP.ts, model = ("ZZZ"))
print(paste("AICc:",round(fit2$aicc,2), "BIC:", round(fit2$bic,2)))

## [1] "AICc: 3579.14 BIC: 3597.26"

#Arima has performed better based on the Information Criterion

#Validation testing
splitDate <- GNP$observation_date[nrow(GNP)] - months(6)
GNP.train <- GNP[GNP$observation_date <= splitDate,]
GNP.valid <- GNP[GNP$observation_date >= splitDate,]

GNP.ts.train <- ts(GNP.train$GNP, start = c(1947,1), end = c(2018,2), frequency = 4)
GNP.ts.valid <- ts(GNP.valid$GNP, start = c(2018,2), end = c(2019,1), frequency = 4)

fit3 <- auto.arima(GNP.ts.train)
fit4 <- ets(GNP.ts.train, model = 'ZZZ')

pred3 <- forecast(fit3, h = 4)
pred4 <- forecast(fit4, h = 4)

print(paste("ARIMA RMSLE:",round(rmsle(pred3$mean, GNP.ts.valid),4)))

## [1] "ARIMA RMSLE: 0.0132"

print(paste("ETS RMSLE:",round(rmsle(pred4$mean, GNP.ts.valid),4)))

## [1] "ETS RMSLE: 0.0179"

#ARIMA has also performed better in our validation method.