ArarForecast

Time series forecasting using ARAR Algorithm. Ref: Introduction to Time Series and Forecasting, Chapter: 10.1.The ARAR Algorithm (Peter J. Brockwell Richard A. Davis (2016) )

Installation

using Pkg
 # dev version
Pkg.add(url = "https://github.com/Akai01/ArarForecast.jl.git")

Usage

Load packages

using CSV
using Downloads
using DataFrames
using TimeSeries
using Dates
using ArarForecast

Download and load the data

dta = CSV.File(Downloads.download("https://raw.githubusercontent.com/Akai01/example-time-series-datasets/main/Data/AirPassengers.csv")) |> DataFrame;
train = filter(row -> row["ds"] < Date(1960,1,1), dta);
test = filter(row -> row["ds"] >= Date(1960,1,1), dta);

Create a TimeArray:

train = (date = train[:,"ds"], data = train[:, "y"]);
train = TimeArray(train; timestamp = :date);
test = (date = test[:,"ds"], test = test[:, "y"]);
test = TimeArray(test; timestamp = :date);
length(test)

## 12

There are different ways to create a TimeArray see TimeSeries.jl package.

Forecasting

fc = arar(;y = train, h = 12, freq = Month, level = [80, 95]);
typeof(fc)

## ArarForecast.Forecast

Plot the Forecast Object

p = ArarForecast.plot(;object = fc)

## Plot{Plots.GRBackend() n=6}

using Plots
Plots.plot(p, test)

The accuracy

accuracy(fc, test)

## (me = [-11.275515996902792], mae = [13.065720800742184], mape = [2.8577953269438026], mdae = [8.718991617233343], rmse = 18.21764747304158)

Benchmark with R forecast package auto.arima

Load the data in and create a ts object

library(magrittr)
dta = read.csv("https://raw.githubusercontent.com/Akai01/example-time-series-datasets/main/Data/AirPassengers.csv")%>%
  dplyr::mutate(ds = as.Date(ds))
head(dta)

##           ds   y
## 1 1949-01-31 112
## 2 1949-02-28 118
## 3 1949-03-31 132
## 4 1949-04-30 129
## 5 1949-05-31 121
## 6 1949-06-30 135

train <- dta%>%dplyr::filter(ds < as.Date("1960-01-01"))

train_ts <- train%>%dplyr::select(-ds)%>%
  ts(start = c(1949, 1), frequency = 12)

test <- dta%>%dplyr::filter(ds >= as.Date("1960-01-01"))

test_ts <- test%>%dplyr::select(-ds)%>%
  ts(start = c(1960, 1), frequency = 12)

Train and forecast 12 months ahead:

fc <- forecast::auto.arima(train_ts)%>%
  forecast::forecast(h = 12)

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

Plot the forecast

forecast::autoplot(fc) + forecast::autolayer(test_ts)

forecast::accuracy(fc$mean, test_ts)

##                 ME    RMSE      MAE       MPE     MAPE       ACF1 Theil's U
## Test set -16.98639 23.9317 18.52768 -3.933491 4.182395 0.04802038 0.5336134

Accuracy Benchmark: R forecast::auto.arima 4.18 vs Julia ArarForecast 2.85

How does the ARAR algorithm Work?

Memory Shortening

The ARAR algorithm applies a memory-shortening transformation if the underlying process of a given time series Y_t, t = 1, 2, ..., n is “long-memory” then it fits an autoregressive model.

The algorithm follows five steps to classify Y_t and take one of the following three actions:

• L: declare Y_t as long memory and form Y_t by Ỹ_t = Y_t − ϕ̂Y_{t − τ̂}
• M: declare Y_t as moderately long memory and form Y_t by Ỹ_t = Y_t − ϕ̂₁Y_t − 1 − ϕ̂₂Y_t − 2
• S: declare Y_t as short memory.

If Y_t declared to be L or M then the series Y_t is transformed again until. The transformation process continuous until the transformed series is classified as short memory. However, the maximum number of transformation process is three, it is very rare a time series require more than 2.

• 1. For each τ = 1, 2, ..., 15, we find the value ϕ̂(τ̂ ) of ϕ̂ that minimizes then define Err(τ) = ERR(ϕ(τ̂ ), τ) and choose the lag τ̂ to be the value of τ that minimizes Err(τ).
• 2. If Err(τ̂) ≤ 8/n, Y_t is a long-memory series.
• 3. If ϕ̂(τ̂) ≥ 0.93 and τ̂ > 2, Y_t is a long-memory series.
• 4. If ϕ̂(τ̂) ≥ 0.93 and τ̂ = 1 or 2, Y_t is a long-memory series.
• 5. If ϕ̂(τ̂) < 0.93, Y_t is a short-memory series.

Subset Autoregressive Model:

In the following we will describe how ARAR algorithm fits an autoregressive process to the mean-corrected series X_t = S_t − S̄, t = k + 1, ..., n where S_t, t = k + 1, ..., n is the memory-shortened version of Y_t which derived from the five steps we described above and S̄ is the sample mean of S_k + 1, ..., S_n.

The fitted model has the following form:

X_t = ϕ₁X_t − 1 + ϕ₁X_t − l1 + ϕ₁X_t − l1 + ϕ₁X_t − l1 + Z

where Z ∼ WN(0,σ²). The coefficients ϕ_j and white noise variance σ² can be derived from the Yule-Walker equations for given lags l₁, l₂, and l₃:

and σ² = γ̂(0)[1−ϕ₁ρ̂(1)] − ϕ_l1ρ̂(l₁)] − ϕ_l2ρ̂(l₂)] − ϕ_l3ρ̂(l₃)], where γ̂(j) and ρ̂(j), j = 0, 1, 2, ..., are the sample autocovariances and autocorelations of the series X_t.

The algorithm computes the coefficients of ϕ(j) for each set of lags where 1 < l₁ < l₂ < l₃ ≤ m where m chosen to be 13 or 26. The algorithm selects the model that the Yule-Walker estimate of σ² is minimal.

Forecasting

If short-memory filter found in first step it has coefficients Ψ₀, Ψ₁, ..., Ψ_k(k≥0) where Ψ₀ = 1. In this case the transforemed series can be expressed as

where Ψ(B) = 1 + Ψ₁B + ... + Ψ_kB^k is polynomial in the back-shift operator.

If the coefficients of the subset autoregression found in the second step it has coefficients ϕ₁, ϕ_l1, ϕ_l2 and ϕ_l3 then the subset AR model for X_t = S_t − S̄ is

where Z_t is a white-noise series with zero mean and constant variance and ϕ(B) = 1 − ϕ₁B − ϕ_l1B^l₁ − ϕ_l2B^l₂ − ϕ_l3B^l₃. From equation (1) and (2) one can obtain

where ξ(B) = Ψ(B)ϕ(B).

Assuming the fitted model in equation (3) is an appropriate model, and Z_t is uncorrelated with Y_j, j < t ∀t ∈ T, one can determine minimum mean squared error linear predictors P_nY_n + h of Y_n + h in terms of 1, Y₁, ..., Y_n for n > k + l₃, from recursions

with the initial conditions P_nY_n + h = Y_n + h, for h ≤ 0.

Ref: Brockwell, Peter J, and Richard A. Davis. Introduction to Time Series and Forecasting. Springer (2016)