Getting started with msarima
In msarima: Multi-Seasonal ARIMA Models

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$.

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.

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.

sacf(MAD_m40, lag.max = 200)

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

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:

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:

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:

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:

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:

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:

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:

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:

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:

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. https://doi.org/10.1002/for.3283

Grigoletto, M. and Lisi, F. (2026). Portmanteau goodness-of-fit tests in the presence of simple, mixed, and multiple seasonal autocorrelation. arXiv. https://doi.org/10.48550/arXiv.2609.32381