README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r setOptions, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# arima2

<!-- badges: start -->
[![R-CMD-check](https://github.com/jeswheel/arima2/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/jeswheel/arima2/actions/workflows/R-CMD-check.yaml)
[![CRAN status](https://www.r-pkg.org/badges/version/arima2)](https://CRAN.R-project.org/package=arima2)
<!-- badges: end -->

The goal of `arima2` is to provide a set of tools to aid in the analysis of time series data in `R`. 
One such function is `arima2::arima`, which provides an interface to estimating Auto Regressive Integrated Moving Average (ARIMA) models using a random-restart algorithm. 
This function improves on the functionality of the `stats::arima` method, as it has the potential to increase the likelihood of the final output model.
By design, the function cannot result in models with lower likelihoods than that of the `stats::arima` function. 
The potential for increased model likelihoods is obtained at the cost of computational efficiency. 
The function is approximately $n$ times slower than the `stats::arima` function, where $n$ is the number of random restarts.
The benefit of trying multiple restarts becomes smaller as the number of available observations increases. 
Because the estimation of ARIMA models takes only a fraction of a second on relatively small data sets (less than a thousand observations), we are of the opinion that potential to increase model likelihoods is well worth this computational cost.
The `arima` function implementation relies heavily on the source code of the `stats::arima` function. 

## Installation

``` r
# Install from CRAN
install.packages("arima2")
```

You can install the development version of `arima2` from [GitHub](https://github.com/) with:

``` r
# install.packages("devtools")
devtools::install_github("jeswheel/arima2")
```

## Example

This is a basic example which shows you how to solve a common problem:

```{r example}
library(arima2)

set.seed(83131)  

# Get model coefficients from ARMA(2, 2)
coefs <- sample_ARMA_coef(order = c(2, 2))

# Get model intercept 
intercept <- rnorm(1, sd = 50)

# Generate data from ARMA model 
x <- intercept + arima.sim(
  n = 100, 
  model = list(ar = coefs[grepl("^ar[[:digit:]]+", names(coefs))], 
               ma = coefs[grepl("^ma[[:digit:]]+", names(coefs))])
)

# Fit ARMA model using arima2 and stats::arima 
arma2 <- arima(x, order = c(2, 0, 2), max_iters = 20)
arma <- stats::arima(x, order = c(2, 0, 2))
```

In the example above, the resulting log-likelihood of the `stats::arima` function is `r round(arma$loglik, 2)`, and the log-likelihood of the `arima` function is `r round(arma2$loglik, 2)`.
For this particular model and dataset, the random restart algorithm implemented in `arima2` improved the model likelihood by `r round(arma2$loglik - arma$loglik, 2)` log-likelihood units. 

Our package creates a new `S3` object that we call `Arima2`, which extends the `Arima` class of the `stats` package. 
Once the model has been fit, our package includes some features that help diagnose the fitted model using this new child class.
For example, `ARMApolyroots` function will return the AR or MA polynomial roots of the fitted model:

```{r getARMAroots}
ARMApolyroots(arma2, type = 'AR')
ARMApolyroots(arma2, type = 'MA')
```

We have also implemented a `plot.Arima2` function that uses the `ggplot2` package so that we can visualize a fitted model. 
To compare the roots of the model fit using multiple restarts to the model fit using `stats::arima`, I will modify the class of the `arma` object so that it can easily be plotted:

```{r PlotARMAresults}
class(arma) <- c("Arima2", class(arma))

plot(arma)
plot(arma2)
```

Finally, if a user would like help in determining an appropriate number of coefficients, we provide the `aicTable` function.
The package also includes an `aicTable` function, which prints the AIC values for all ARMA$(p, d, q)$, with $p \leq P$, $q \leq Q$, and $d = D$:

```{r getAICtable}
set.seed(443252)
tab_results <- aicTable(x, P = 4, Q = 4, D = 0) 

tab_results |> knitr::kable()

P <- which(tab_results == min(tab_results), arr.ind = TRUE)[1] - 1
Q <- which(tab_results == min(tab_results), arr.ind = TRUE)[2] - 1

print(paste0("p = ", P, "; q = ", Q))
```

For more details about this package, please see our arXiv paper: [arXiv:2310.01198](https://doi.org/10.48550/arXiv.2310.01198)