--- title: "Getting started with msarima" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with msarima} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 8, fig.height = 4 ) ``` Hourly, daily or half-hourly data usually show more than one seasonal cycle at the same time: the number of cars on a motorway repeats itself every 24 hours *and* every 168 hours (a week). `stats::arima()` can accommodate only one of the two. The `msarima` package fits the multiple seasonal model of Lisi and Grigoletto (2025), $$ \phi(B)\,\prod_{j=1}^{k}\Phi_j\!\left(B^{S_j}\right)\, (1-B)^{d}\prod_{j=1}^{k}\left(1-B^{S_j}\right)^{D_j} y_t = \mu + \theta(B)\,\prod_{j=1}^{k}\Theta_j\!\left(B^{S_j}\right) e_t , $$ where $e_t$ is white noise with variance $\sigma^2$ and $\mu$ is a constant, with one autoregressive, differencing and moving average operator for each seasonal period $S_1, \dots, S_k$. ```{r setup, message = FALSE, warning = FALSE} library(msarima) ``` ## The data `MAD_m40` is the hourly number of vehicles recorded by a sensor on the M-40 ring road of Madrid over the eight weeks from 3 February to 29 March 2024. It is an `msts` object with seasonal periods 24 and 168 (1344 observations in all), so those periods do not have to be repeated in every call. ```{r data, fig.alt = "Line chart of vehicles per hour over the first four weeks. A daily cycle of two peaks and a night trough repeats, with lower peaks on the last two days of every week."} head(MAD_m40, 12) plot(window(MAD_m40, end = 4), ylab = "vehicles per hour", main = "M-40, first four weeks (time in weeks)") ``` Both cycles are visible: the daily peak-hour pattern, and the weekly one with the weekends clearly lower. ## Looking at the cycles `sacf()` and `spacf()` are the correlograms of the package. They work like `stats::acf()`, but the lags belonging to each seasonal period are coloured, so the cycles can be read off the plot instead of counted by hand. The periods are taken from the `msts` object. ```{r acf, fig.height = 3.5, fig.alt = "Correlogram up to lag 200. The autocorrelations stay well above the confidence band and peak at the multiples of 24, in orange, and at 168, in blue, decaying very slowly."} sacf(MAD_m40, lag.max = 200) ``` ```{r pacf, fig.height = 3.5, fig.alt = "Partial correlogram up to lag 200. Only the first lags and those around the multiples of 24 and 168 stand out of the confidence band."} spacf(MAD_m40, lag.max = 200) ``` Orange marks the multiples of 24, blue the multiples of 168; `col` changes them. The autocorrelation decays very slowly at both, which suggests one difference at each period. There is a `ggplot2` version of both, `ggsacf()` and `ggspacf()`: ```{r ggacf, fig.height = 3.5, eval = requireNamespace("ggplot2", quietly = TRUE), fig.alt = "The same correlogram in ggplot2, with a legend mapping orange to period 24 and blue to period 168."} ggsacf(MAD_m40, lag.max = 200) ``` ## Fitting the model The non-seasonal part goes in `order`, and each seasonal component is one element of the `seasonal` list, in the same order as the periods: ```{r fit} fit <- msarima(MAD_m40, order = c(1, 0, 2), seasonal = list(c(0, 1, 1), c(0, 1, 1))) fit ``` `summary()` adds standard errors, obtained from the numerical Hessian of the log-likelihood, and the associated tests: ```{r summary} summary(fit) ``` `coef()`, `fitted()`, `residuals()`, `logLik()`, `AIC()` and `BIC()` are all available and behave as usual. When the model contains differencing the series is differenced before estimation, so the first $d + \sum_k D_k S_k$ residuals are `NA`: ```{r nused} c(nobs = fit$nobs, n.used = fit$n.used, leading_NA = sum(is.na(residuals(fit)))) ``` Parameters can be held fixed by supplying `fixed` with one entry per parameter and `NA` where the parameter is to be estimated. The order is AR, seasonal AR, MA, seasonal MA, sigma: ```{r fixed} fit_fixed <- msarima(MAD_m40, order = c(1, 0, 2), seasonal = list(c(0, 1, 1), c(0, 1, 1)), fixed = c(NA, NA, NA, -0.6, NA, NA)) c(free = AIC(fit), fixed = AIC(fit_fixed)) ``` ## Checking the residuals Passing the fitted model to `sacf()` uses the residuals and the seasonal periods of the model, so no argument is needed. What is left at the seasonal lags is what the model has failed to capture: ```{r resid, fig.height = 3.5, fig.alt = "Correlogram of the residuals up to lag 200. The spikes are much smaller than before and mostly inside the confidence band."} sacf(fit, lag.max = 200, main = "Residuals") ``` `msportmanteau()` turns this visual check into a formal test, using the multi-seasonal portmanteau statistics of Grigoletto and Lisi (2026). With the first 8 lags and 5 multiples of each period, 18 lags are tested, leaving 13 degrees of freedom after the 5 estimated coefficients. Neither the Ljung-Box nor the Monti version rejects the absence of residual autocorrelation: ```{r portm} msportmanteau(fit, nlag = 5, lags_short = 8) msportmanteau(fit, type = "Monti", nlag = 5, lags_short = 8) ``` ## Forecasting `forecast()` returns an object of the usual class of the **forecast** package, so it prints and plots like any other. Prediction intervals come from the $\mathrm{MA}(\infty)$ representation of the model, under normal innovations: ```{r forecast, fig.alt = "Line chart of the last two weeks of data followed by 48 hours of forecasts, which repeat the daily shape, inside shaded 80 and 95 per cent prediction bands."} fc <- forecast(fit, h = 48) plot(fc, include = 336, main = "Two days ahead", ylab = "vehicles per hour") ``` ## Simulating `msarima.sim()` generates a series from a model specification, without any data: ```{r sim} set.seed(1) y <- msarima.sim(nsim = 500, order = c(1, 0, 1), seasonal = list(c(0, 0, 1), c(0, 0, 1)), seasonal.period = c(7, 30), ar = 0.6, ma = c(0.4, 0.3, 0.2), sigma2 = 0.5) ``` and `simulate()` does the same from a fitted model, reusing its specification and its estimates: ```{r sim2} y_sim <- simulate(fit, nsim = 168, seed = 1) ``` ## References Lisi, F. and Grigoletto, M. (2025). Multiple seasonal autoregressive integrated moving average models. *Journal of Forecasting*, **44**(6), 2037--2052. Grigoletto, M. and Lisi, F. (2026). Portmanteau goodness-of-fit tests in the presence of simple, mixed, and multiple seasonal autocorrelation. *arXiv*.