Package 'muse'

Initial state values for a fitted pts object.

Description

Returns the smoothed structural states at the first observation, with the "Error" and "Fit" columns dropped. For deterministic components – e.g. the slope under the global (G) trend, where the slope is fixed throughout the horizon – this equals the t = 0 initial value; for stochastic components it is the smoother's estimate at t = 1, a close proxy to the t = 0 initial.

Usage

initials(object, ...)

Arguments

object	A fitted object of class "pts".
...	Unused.

Value

Named numeric vector of initial structural-state values.

---

muse package

Description

Package contains functions implementing Multiple Source of Error state space models for purposes of time series analysis and forecasting.

Details

Package:	muse
Type:	Package
Date:	2024-10-18 - Inf
License:	LGPL-2.1

The following functions are included in the package:

Author(s)

Diego J. Pedregal, [email protected]

Ivan Svetunkov, [email protected]

See Also

adam

Examples

x <- rnorm(100,100,10)

---

pts: Power / Trend / Seasonal state-space model

Description

Estimates a PTS (Power / Trend / Seasonal) state-space model for a univariate time series. This is the user-facing entry point of the muse package and mirrors the calling convention used elsewhere in the smooth family: pts() estimates the model, and forecast.pts produces forecasts from the fitted object without re-estimating.

Usage

pts(data, model = "ZZZ", lags = stats::frequency(data), orders = list(ar
  = 0, ma = 0, select = FALSE), formula = NULL, regressors = c("use"),
  outliers = c("ignore", "use", "select"), level = 0.99, ic = c("AICc",
  "BICc", "BIC", "AIC"), lambda_estim = c("likelihood", "guerrero",
  "decomp-guerrero"), biasadj = FALSE, h = 0, holdout = FALSE,
  verbose = FALSE, ...)

Arguments

data	response series. Either a univariate ts / zoo / numeric vector, OR a matrix / data.frame whose first column is the response and whose remaining columns are external regressors (xregs).
model	3-letter PTS specification string. The three positions encode Power / Trend / Seasonal: Power: Z to estimate Box-Cox $\lambda$, or a numeric value (e.g. "0", "0.5", "1"). Trend: Z (auto), N (none / random walk), L (local linear), D (damped / smooth random walk), G (global / deterministic). Seasonal: Z (auto), N (none), D (discrete / linear), T (trigonometric / equal).
lags	seasonal period (default frequency(data)). Scalar for the structural seasonal; may also be a vector c(1, s) mirroring smooth::adam's convention – the last entry is the structural period, the full vector becomes the default lag set for the irregular's ARMA blocks (overridable per-fit via orders$lags).
orders	ARMA / SARMA spec for the irregular component. Three forms: Full list list(ar, ma, select) for a non-seasonal $ARMA(p, q)$ – ar, ma non-negative scalars (default 0); select = TRUE asks the engine to search ARMA orders up to that cap. Numeric shortcut c(p, q) – equivalent to list(ar = p, ma = q, select = FALSE); c(p) is treated as c(p, 0). Seasonal list(ar = c(p, P), ma = c(q, Q), lags = c(1, s)) for $SARMA(p, q)(P, Q)_s$ – ar / ma are length-L vectors paired position-wise with lags, with lags[1] = 1 (non-seasonal block) and lags[2] = s (seasonal block). If orders$lags is omitted the default falls back to the top-level lags argument (or c(1, frequency(data))). The seasonal SARMA polynomial is multiplied internally, so the BFGS only optimises the free $\phi_i, \Phi_j, \theta_i, \Theta_j$ coefficients. select = TRUE runs a grid search over every $(p', q', P', Q')$ tuple with $0 \le p' \le$ ar[1] and so on, and picks the candidate with the lowest ic. PTS has no differencing – any orders$i supplied is silently ignored.
formula	optional formula response ~ x1 + x2 + ...; only meaningful when data is a matrix or data.frame. Used to pick the response column + xreg columns explicitly.
regressors	handling of xregs. Currently only "use" (apply all supplied xregs as fixed-coefficient covariates). Adam's "select" and "adapt" modes are not yet implemented.
outliers	what to do about outliers, mirroring adam()'s interface: "ignore" (default – fit ignores the possibility of outliers) or "use" (run the engine's outlier detector once after the structural fit, classify each event as AO / LS / SC, and refit with the detected events as fixed regressor dummies). Adam's "select" mode (IC-pruning of detected dummies) is not yet supported – passing it errors with a clear message. When outliers = "use" the detected events are returned on the fitted object as $outliersDetected (a data frame with time and type columns) and the corresponding dummy coefficients appear in coef(m) as Beta(...) entries.
level	confidence level driving the outlier z-score threshold (default 0.99). Translated to a two-sided z via qnorm((1 + level) / 2): 0.99 -> ~= 2.576, 0.95 -> ~= 1.960. The z drives the AO threshold; LS / SC scale proportionally to preserve the engine's relative stiffness (LS ~= 1.087*z, SC ~= 1.304*z). Ignored when outliers = "ignore".
ic	information criterion for automatic model selection; one of "AICc" (default), "BICc", "BIC", "AIC". Matches the adam option set; AICc is the default, as in adam.
lambda_estim	how the Box-Cox power $\lambda$ is chosen when the power slot of model is "Z"; one of: "likelihood" (default) – estimate $\lambda$ jointly with the structural parameters by maximum likelihood in the engine. "guerrero" – the classical Guerrero (1993) variance- stabilisation screen on raw season-length blocks. "decomp-guerrero" – Guerrero on an msdecompose-smoothed trend (the former default). Ignored when a numeric power is supplied (e.g. "0.5ZZ").
biasadj	logical (default FALSE). Point forecasts are the back-transformed conditional median $g^{-1}(\mu)$. When TRUE, a second-order bias correction is applied so the point forecast approximates the conditional mean (as in forecast::InvBoxCox(biasadj = TRUE)). Prediction-interval quantiles are unaffected. Has no effect at $\lambda = 1$.
h	forecast horizon. If h > 0 a forecast is computed at fit time and cached on the object; forecast(object, h) can later recompute for a different horizon cheaply.
holdout	logical. If TRUE and h > 0, the last h observations of data are withheld from estimation and returned in $holdout for later accuracy assessment.
verbose	logical: print intermediate optimisation output.
...	advanced / undocumented passthroughs. Supported keys: B - numeric vector of starting values for the optimiser (natural-scale variances, in the order returned in $B by a default fit). Mirrors the same hatch in smooth::adam(). The optimised vector is returned in the $B slot of the output regardless of whether the user supplied one.

Value

An object of class c("pts", "smooth"). Slot names mirror smooth::adam()'s return list where the concept is shared; pts-only extensions are flagged below.

• Inputs / spec: y, model, modelUC*, lags, lambda*

• Parameters: B (estimated parameter vector), covp* (parameter covariance), nParam – an adam-style 2 x 5 matrix (rows Estimated / Provided; columns nParamInternal, nParamXreg, nParamOccurrence, nParamScale, nParamAll). nparam() returns the [Estimated, nParamAll] cell. Estimated structural initials (level/slope, cycle, seasonal states) are folded into nParamInternal, exactly as smooth::adam does

• In-sample fit: fitted, residuals, states plus pts-specific comp* (additive BC-scale decomposition with Error/Fit columns)

• Cached forecast: forecast (original scale, if h > 0) and forecast_args* for cheap re-forecasting

• Likelihood + diagnostics: logLik, table* (C++ validation text)

• Scalars read by plot.smooth / diagnostics: distribution = "dnorm", loss = "likelihood", occurrence = NULL, holdout

• Bookkeeping: call, timeElapsed

AIC / AICc / BIC / BICc are derived on demand via the methods, not stored on the object. (* = pts-specific extension.)

Author(s)

Diego J. Pedregal, [email protected]

Ivan Svetunkov, [email protected]

See Also

forecast.pts

Examples

# Automatic model selection (Power / Trend / Seasonal) on monthly data
model <- pts(AirPassengers, model = "ZZZ", h = 12, holdout = TRUE)
model

# A fixed specification: no transform, local-linear trend, trigonometric
# seasonality, with a 12-step forecast
fixedModel <- pts(AirPassengers, model = "1LT", h = 12)
forecast(fixedModel, h = 12)

---

S3 methods for objects of class pts

Description

Standard accessors and printing for fitted pts objects. forecast(object, h) generates forecasts from the fitted parameters without re-running the optimiser; predict(object) returns the in-sample fitted values.

print.pts mirrors smooth::print.adam (smooth/R/adam.R:5862) section by section, substituting the PTS-specific MSOE concepts where ETS-specific ones do not apply: the "initialisation" line becomes the Box-Cox $\lambda$ block, and the "Persistence vector g" block becomes the MSOE innovation variances (the same variances are otherwise the structural pieces of the B parameter vector). The C++ validation diagnostics live on $cppOutput should the user want them.

forecast.pts uses the C++ forecastOnly entry point: it skips re-estimation and just propagates the Kalman filter h steps forward from the fitted state, so changing h is cheap. interval selects the variance source:

"prediction" (default)

state propagation + future shocks (the engine's yForV). Bands the next observation.

"confidence"

conditional-mean variance. In a state-space model $var(E[y_{t+h}\,|\,obs]) = var(y_{t+h}\,|\,obs) - \sigma^2$, where $\sigma^2$ is the BC-scale residual variance: we read it straight off the fitted scale, no reforecast needed.

"simulated"

empirical quantiles of nsim forward paths from simulate.pts(). Returns the path matrix in $scenarios when scenarios = TRUE.

"none"

no bands; lower = upper = mean.

level accepts a vector; lower / upper are then $h \times nLevels$ matrices (or $1 \times nLevels$ when cumulative = TRUE). side produces two-sided / upper-only / lower-only intervals – the absent side is set to $\mp\infty$ on the BC scale and back-transformed to the original-scale support boundary (0 for $\lambda > 0$, $-\infty$ for the identity transform). cumulative = TRUE collapses the horizon into one total – exact for interval = "simulated" (sum within each path); the engine does not expose cross-step state covariance, so the other intervals fall back to simulation totals.

Usage

## S3 method for class 'pts'
print(x, digits = 4, ...)

## S3 method for class 'pts'
plot(x, which = c(1, 2, 4, 6), ...)

## S3 method for class 'pts'
fitted(object, ...)

## S3 method for class 'pts'
residuals(object, ...)

## S3 method for class 'pts'
coef(object, ...)

## S3 method for class 'pts'
vcov(object, ...)

## S3 method for class 'pts'
nobs(object, all = FALSE, ...)

## S3 method for class 'pts'
logLik(object, ...)

## S3 method for class 'pts'
predict(object, newdata = NULL, ...)

## S3 method for class 'pts'
forecast(object, h = 10, newdata = NULL,
  interval = c("prediction", "confidence", "simulated", "none"),
  level = 0.95, side = c("both", "upper", "lower"), cumulative = FALSE,
  nsim = NULL, scenarios = FALSE, biasadj = NULL, ...)

## S3 method for class 'pts'
initials(object, ...)

Arguments

x, object	A fitted object of class "pts".
digits	number of significant digits used in printed output.
...	further arguments passed to underlying generics.
which	integer vector of plot panels to draw; passed through to plot.smooth (see ?smooth::plot.adam). Defaults to c(1, 2, 4, 6).
all	logical; if TRUE the holdout sample is included in the observation count. Default FALSE.
newdata	unused (reserved for forecast paths with new regressors).
h	forecast horizon (steps ahead).
interval	interval type; one of "prediction", "confidence", "simulated", or "none".
level	confidence level for prediction intervals (default 0.95).
side	one of "both", "upper", "lower" – selects two-sided vs one-sided intervals.
cumulative	if TRUE, return the cumulative h-step total.
nsim	number of simulated paths when interval = "simulated" or under the cumulative fallback (default 10000).
scenarios	if TRUE and interval = "simulated", return the simulated path matrix as $scenarios.
biasadj	point-forecast back-transform: FALSE (median, the default) or TRUE (bias-corrected mean). Defaults to the value the model was fitted with (object$biasadj).

Value

The methods return the following:

• fitted, predict – the in-sample fitted values as a ts / numeric vector on the original (back-transformed) scale.

• residuals – the model residuals as a ts / numeric vector on the Box-Cox scale.

• coef – a named numeric vector of the estimated parameters (structural variances plus any ARMA / damping coefficients).

• vcov – the parameter variance-covariance matrix.

• nobs – an integer, the number of in-sample observations (plus the holdout when all = TRUE).

• nparam – an integer, the number of estimated parameters (the degrees of freedom used by the information criteria).

• logLik – an object of class "logLik": the maximised log-likelihood, with df and nobs attributes.

• sigma, extractSigma – a numeric scalar, the residual standard deviation; extractScale – the maximum-likelihood scale of the assumed distribution.

• actuals – the original response series.

• modelType – a character string with the PTS model code (e.g. "PTS(0,N,T)"); modelName – a human-readable description; errorType – the character error type ("A", additive on the Box-Cox scale).

• lags – the integer seasonal period; orders – a list of the irregular ARMA orders (ar, ma, lags); initials – a named numeric vector of the estimated initial structural-state values.

• forecast – an object of class "pts.forecast": a list with the point forecast $mean and the interval bounds $lower / $upper (original scale), together with the forecast variance, level, and interval metadata.

• print and plot are called for their side effects (console output and a diagnostic plot, respectively) and invisibly return their first argument.

Author(s)

Diego J. Pedregal, [email protected]

Ivan Svetunkov, [email protected]

References

• Granger, C. W. J., & Newbold, P. (1976). Forecasting transformed series. Journal of the Royal Statistical Society: Series B (Methodological), 38(2), 189-203. doi:10.1111/j.2517-6161.1976.tb01585.x

• Pankratz, A., & Dudley, U. (1987). Forecasts of power-transformed series. Journal of Forecasting, 6(4), 239-248. doi:10.1002/for.3980060403

• Guerrero, V. M. (1993). Time-series analysis supported by power transformations. Journal of Forecasting, 12(1), 37-48. doi:10.1002/for.3980120104

Examples

model <- pts(AirPassengers, model = "1LT", h = 12, holdout = TRUE)
print(model)
fitted(model)
residuals(model)
# forecast 12 steps ahead with 95% prediction intervals
forecast(model, h = 12, interval = "prediction", level = 0.95)

---