Neural Network modeling for time series forecast of US Trade Balance (2016 - 2020). Data obtained from U.S. Census Bureau.

Load Liabilities

library(fpp3)
## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──
## ✔ tibble      3.2.1     ✔ tsibble     1.1.3
## ✔ dplyr       1.1.2     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.0     ✔ feasts      0.3.1
## ✔ lubridate   1.9.2     ✔ fable       0.3.3
## ✔ ggplot2     3.4.2     ✔ fabletools  0.3.4
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
library(readxl)

Get data

mydata <- read_excel('/Users/adrianjones/Documents/Forecasting/Week6/Discussion6/Trade.xlsx')

library(zoo)
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
mydata$Period <- as.yearmon(mydata$Period, format = "%b-%Y")

mydata <- mydata[133:192,] #Years 2016 - 2020

Time series + train/test

myts <- mydata %>%
  mutate(Period = yearmonth(Period)) %>%
  as_tsibble(index = Period)

ggplot(myts, aes(x=Period, y=Value)) +
  geom_line(aes(y=Value)) +
  ggtitle("US Trade Balance")

# Train and Test sets
train <- myts[1:48,]
test <- myts[49:nrow(myts),]

Time Series Plots

myts%>%gg_season(Value)

myts%>%gg_subseries(Value)

myts%>%gg_lag(Value)

myts%>%ACF(Value)%>%autoplot()

Transformation

lambda <- myts |> features(Value, features = guerrero) |> pull(lambda_guerrero)
## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value

## Warning in optimise(lambda_coef_var, c(lower, upper), x = x, .period =
## max(.period, : NA/Inf replaced by maximum positive value
myts |> autoplot(box_cox(Value, lambda)) +
  labs(y = "",title = latex2exp::TeX(paste0(
    "Transformed US Trade Balance $\\lambda$ = ", round(lambda,2))))

# Decomposition
comp=myts%>%model(stl=STL(Value))%>%components()
comp|>as_tsibble() |>autoplot(Value)+  geom_line(aes(y=trend), colour = "red")+
  geom_line(aes(y=season_adjust), colour = "blue")

comp%>%autoplot()

Models (This takes a bit to run)

m1 <- train |> model(Seasonal_model = SNAIVE(Value))
m2 <- train |> model(Neural_model = NNETAR((Value) ~ trend() + season()))
m3 <- train |> model(ETS_model = ETS(Value))
m4 <- train |> model(Ensemble = (SNAIVE(Value)) + NNETAR((Value)~trend()+season()) + ETS(Value)/3)

Forecast reports

for (i in 1:4){report(eval(parse(text=paste0('m',i))))}
## Series: Value 
## Model: SNAIVE 
## 
## sigma^2: 29066840.3071 
## Series: Value 
## Model: NNAR(2,1,8)[12] 
## Transformation: (Value) 
## 
## Average of 20 networks, each of which is
## a 15-8-1 network with 137 weights
## options were - linear output units 
## 
## sigma^2 estimated as 1.028
## Series: Value 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.5632494 
## 
##   Initial states:
##      l[0]
##  -39847.4
## 
##   sigma^2:  7788298
## 
##      AIC     AICc      BIC 
## 951.4452 951.9906 957.0588 
## Series: Value 
## Model: COMBINATION 
## Combination: Value + Value
## 
## ==========================
## 
## Series: Value 
## Model: COMBINATION 
## Combination: Value + Value
## 
## ==========================
## 
## Series: Value 
## Model: SNAIVE 
## 
## sigma^2: 29066840.3071 
## 
## Series: Value 
## Model: NNAR(2,1,8)[12] 
## Transformation: (Value) 
## 
## Average of 20 networks, each of which is
## a 15-8-1 network with 137 weights
## options were - linear output units 
## 
## sigma^2 estimated as 1.477
## 
## 
## Series: Value 
## Model: COMBINATION 
## Combination: Value * 0.333333333333333
## 
## ======================================
## 
## Series: Value 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.5632494 
## 
##   Initial states:
##      l[0]
##  -39847.4
## 
##   sigma^2:  7788298
## 
##      AIC     AICc      BIC 
## 951.4452 951.9906 957.0588

Plots

myplot=function(model){
  model|>forecast(test)|>autoplot(myts)+geom_line(aes(y=.fitted, col=.model),data=augment(model))+ggtitle(names(model))}

for (i in 1:4){print(myplot(eval(parse(text=paste0('m',i)))))}
## Warning: Removed 12 rows containing missing values (`geom_line()`).

## Warning: Removed 12 rows containing missing values (`geom_line()`).

## Warning: Removed 12 rows containing missing values (`geom_line()`).

Accuracy Metrics

tmp3=c(rep(0,10))
mymat=matrix(rep(0, 4*nrow(train)), ncol=4)

for (i in 1:4){
  tmp=eval(parse(text=paste0('m',i)))%>%forecast(test)
  tmp2=tmp%>%accuracy(test,measures = list(point_accuracy_measures))
  tmp3=rbind(tmp2, tmp3)
}

tmp3=tmp3[tmp3$RMSE>0,]
tmp3$MASE=tmp3$RMSSE=tmp3$.type=NULL
tmp3%>%arrange(RMSE)
## # A tibble: 4 × 7
##   .model              ME   RMSE    MAE    MPE  MAPE  ACF1
##   <chr>            <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>
## 1 Neural_model    -6375.  8291.  7256.   10.5  12.7 0.514
## 2 Seasonal_model  -7790. 12810. 10515.   11.8  18.2 0.758
## 3 ETS_model      -11631. 14046. 11991.   19.6  20.4 0.692
## 4 Ensemble        55334. 55777. 55334. -106.  106.  0.675