Introduction

The muse package fits PTS (Power / Trend / Seasonal) state-space models — single-source-of-error structural time-series models whose components mirror the ETS/ADAM taxonomy. All estimation happens through a single user-facing function, pts(), which delegates to a C++ engine (Rcpp + RcppArmadillo) running a Kalman filter / smoother under a concentrated likelihood.

pts() returns an object of class c("pts", "smooth"), so all generics shared with the smooth / greybox ecosystem (forecast(), plot(), accuracy(), AIC(), BIC(), simulate(), …) work out of the box.

This vignette walks through the main features on two datasets:

• AirPassengers — the classic monthly airline-passenger series, useful for showing the Box-Cox power transformation and automatic model selection.
• Seatbelts — UK road-casualty data including a structural break (the 1983 seat-belt law), useful for showing external regressors and engine-side outlier detection.

library(muse)
#> Loading required package: greybox
#> Package "greybox", v2.0.9.41002 loaded.
#> Loading required package: smooth
#> This is package "smooth", v4.5.0

Quick start: AirPassengers

The simplest call lets pts() choose everything automatically:

air <- pts(AirPassengers, model = "ZZZ", h = 12, holdout = TRUE)
air
#> Time elapsed: 0.42 seconds
#> Model estimated using pts() function: PTS(2,L,T)
#> With Box-Cox lambda = 2
#> Periods: 12.0 / 6.0 / 4.0 / 3.0 / 2.4 / 2.0
#> Distribution assumed in the model: Normal
#> Loss function type: likelihood; Loss function value: 483.11
#> Variance parameters:
#>               Variance Proportion
#> Level          25270.0   0.087220
#> Slope          68160.0   0.235300
#> Seas(All) (*) 195400.0   0.674400
#> Irregular        889.8   0.003071
#> (*) concentrated out
#> 
#> Sample size: 132
#> Number of estimated parameters: 17
#> Number of degrees of freedom: 115
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 1000.220 1005.588 1049.228 1062.334

The model string is a compact three-character spec encoding Power / Trend / Seasonal in that order; the three Zs ask the engine to choose:

• the Box-Cox power λ (jointly with the state-space parameters),
• the trend type,
• the seasonal type.

Because we asked for holdout = TRUE with h = 12, the last 12 months were withheld from estimation and stashed on the object for later accuracy evaluation.

A look at the in-sample fit and residuals:

plot(air, which = c(1, 2))

To pull the structural decomposition (level / slope / seasonal / irregular):

plot(air, which = 12)

A forecast over the holdout window:

fc <- forecast(air, h = 12)
plot(fc)

…and accuracy against the held-out values:

accuracy(air)
#>            ME           MAE           MSE           MPE          MAPE 
#>  -3.141551008  15.624135272 416.959651342  -0.011783682   0.033972940 
#>           sCE          sMAE          sMSE          MASE         RMSSE 
#>  -0.143617905   0.059522233   0.006051464   0.648735886   0.651714094 
#>          SAME          rMAE         rRMSE          rAME     asymmetry 
#>   0.130441579   0.205580727   0.198293621   0.044143574  -0.119834000 
#>          sPIS 
#>   0.873553552

The PTS model string

The model argument is a compact three-character string of the form "PTS" (Power-Trend-Seasonal): position 1 is the Box-Cox λ, position 2 the trend letter, position 3 the seasonal letter. Numeric λ values can take more than one character ("0.5LT", "-0.3LD") — the trend and seasonal letters are always read from the last two positions. The fitted object’s $model field prints in a more readable form (e.g. "PTS(-0.3, G, D)") but that pretty-printed string is not accepted as input.

air_log <- pts(AirPassengers, model = "0LT", h = 12, holdout = TRUE)
air_log$lambda
#> [1] 0

Automatic selection

Setting any letter to Z triggers an internal search over the admissible candidates. The information criterion is set via ic (default "AICc"):

pts(AirPassengers, model = "ZZZ", h = 12, ic = "BIC")

ARMA on the irregular component

The irregular component can itself be an ARMA / SARMA process, controlled by orders:

# ARMA(2, 1) on the irregular
pts(AirPassengers, model = "ZZZ", orders = c(2, 1), h = 12)

# Seasonal SARMA(1, 1)(1, 1)_12 — non-seasonal + seasonal blocks
pts(AirPassengers, model = "ZZZ",
    orders = list(ar = c(1, 1), ma = c(1, 1), lags = c(1, 12)), h = 12)

# Automatic ARMA search up to the supplied caps
pts(AirPassengers, model = "ZZZ",
    orders = list(ar = 2, ma = 2, select = TRUE), h = 12)

When orders$select = TRUE muse runs a nested PTS-outer / ARMA-inner selection: for every structural candidate it scores the best ARMA order found on the residuals, and picks the joint (structure, ARMA) pair with the lowest IC.

Forecasting

forecast() is the workhorse for producing point forecasts and prediction intervals. The mean is always returned on the original scale (back-transformed from Box-Cox); intervals are built by endpoint-transforming the BC-scale bounds.

fc <- forecast(air, h = 12, interval = "prediction", level = c(0.80, 0.95))
plot(fc)

Three interval types are supported:

cumulative = TRUE collapses the forecast horizon into a single total — exact for "simulated", an approximation otherwise.

Outlier detection

pts() can detect and absorb outliers in a single call, mirroring the adam-style outliers / level API:

y <- AirPassengers
y[100] <- 3 * y[100]            # inject an obvious additive spike

m_out <- pts(y, model = "ZZZ", outliers = "use", level = 0.99)
m_out$outliersDetected
#>   time type
#> 1   23   AO
#> 2  100   AO

The detected events are classified as:

• AO — additive outlier (one-off spike),
• LS — level shift (step change),
• SC — slope change.

level controls the underlying z-score threshold via qnorm((1 + level) / 2) — at 0.99 ≈ 2.576. Detected outliers are injected as fixed dummy regressors (AO<t>, LS<t>, SC<t>); their coefficients show up in coef(m_out) and print() reports a one-line summary of how many of each type were absorbed.

Setting outliers = "ignore" (the default) skips the detection step entirely.

External regressors

pts() accepts a data.frame / matrix / mts whose first column is the response and whose remaining columns are external regressors. On the Seatbelts dataset we can use kms, PetrolPrice, and law (the 1983 seat-belt law indicator) as regressors for monthly drivers (drivers killed or seriously injured):

sb <- Seatbelts[, c("drivers", "kms", "PetrolPrice", "law")]
m_sb <- pts(sb, model = "ZZZ", h = 12, holdout = TRUE)
m_sb
#> Time elapsed: 0.81 seconds
#> Model estimated using pts() function: PTS(0.5,N,D)
#> With Box-Cox lambda = 0.5
#> Distribution assumed in the model: Normal
#> Loss function type: likelihood; Loss function value: 1144.596
#> Variance parameters:
#>           Variance Proportion
#> Level        1.243     0.2260
#> Seas (*)     1.243     0.2260
#> Irregular    3.012     0.5479
#> (*) concentrated out
#> 
#> Regressor coefficients:
#>    Beta(1)    Beta(2)    Beta(3) 
#>  2.163e-04 -1.190e+02 -9.816e+00 
#> 
#> Sample size: 180
#> Number of estimated parameters: 21
#> Number of degrees of freedom: 159
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 2331.191 2337.040 2398.244 2413.428

For point forecasting with regressors you must supply newdata covering the forecast horizon:

# Holdout values of the regressors are stashed on $holdout for accuracy
fc_sb <- forecast(m_sb, h = 12, newdata = tail(Seatbelts, 12))
plot(fc_sb)

A pure structural-break example

The law column makes Seatbelts a natural showcase for engine-side outlier detection: the 1983 intervention should look like a level shift to a structural model that does not have it as a regressor.

m_sb_out <- pts(Seatbelts[, "drivers"], model = "ZZZ",
                outliers = "use", level = 0.99)
m_sb_out$outliersDetected
#> [1] time type
#> <0 rows> (or 0-length row.names)

The engine is expected to flag a level shift around February 1983 (observation 170 of 192).

Holdout evaluation and accuracy

The combination of holdout = TRUE plus h > 0 reserves the tail of the series for evaluation. accuracy() then refits-free scores the forecast against the holdout:

accuracy(air)
#>            ME           MAE           MSE           MPE          MAPE 
#>  -3.141551008  15.624135272 416.959651342  -0.011783682   0.033972940 
#>           sCE          sMAE          sMSE          MASE         RMSSE 
#>  -0.143617905   0.059522233   0.006051464   0.648735886   0.651714094 
#>          SAME          rMAE         rRMSE          rAME     asymmetry 
#>   0.130441579   0.205580727   0.198293621   0.044143574  -0.119834000 
#>          sPIS 
#>   0.873553552

The metrics come from greybox::measures() so the full suite (MAE, MASE, RMSE, MAPE, sCE, …) is available.

Simulation

simulate() produces fan trajectories from the fitted model:

set.seed(42)
sim <- simulate(air, nsim = 200, h = 12)
str(sim, max.level = 1)
#> List of 5
#>  $ data : Time-Series [1:132, 1:200] from 1949 to 1960: 112 125 125 149 107 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ model: chr "PTS(2,L,T)"
#>  $ nsim : int 200
#>  $ obs  : int 132
#>  $ seed : NULL
#>  - attr(*, "class")= chr [1:2] "pts.sim" "smooth.sim"

In-sample replay starts from the initial state estimated by the filter; out-of-sample paths propagate forward with state and observation disturbances drawn from the fitted noise distribution.

Inspecting the fit

A range of standard accessors works on pts objects:

summary(air)         # Coefficient table + variance proportions
coef(air)            # Estimated parameter vector
vcov(air)            # Parameter covariance matrix
confint(air)         # Wald confidence intervals
logLik(air); AIC(air); BIC(air)
fitted(air); residuals(air)
rstandard(air); rstudent(air)
nparam(air); nobs(air); sigma(air)
modelType(air); orders(air); lags(air); errorType(air)

Further reading

• ?pts — full argument documentation including the orders shortcuts and the outliers / level interface.
• ?forecast.pts — interval types, side, cumulative, and the scenarios argument.
• The ARCHITECTURE.md file in the package root for a tour of the R ↔︎ C++ boundary and the internal state-space machinery.