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()

Model

fit <- train %>% 
  model(NNETAR((Value) ~ trend() + season()))
fit %>%
  forecast(h=12) %>%
  autoplot(myts) +
  labs(x= "Period", y= "Trade Balance", title = "Monthly US Trade Balance")

Forecast report

forecast_results <- fit %>%
  forecast(h=12)

# Report
report(fit)
## 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.606

Accuracy Metrics

accuracy(forecast_results, test) #measure accuracy of forecast against held out sample
## # A tibble: 1 Ă— 10
##   .model                  .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>                   <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NNETAR((Value) ~ trend… Test  -5807. 7845. 6912.  9.47  12.1   NaN   NaN 0.559