Time series analysis – method notes #

The decisions behind the time series analysis module – which estimator, threshold, refusal or fallback stands behind each card, and the numbers that validate them, each computed on the airline passenger table printed on the time series concepts page. Sections follow the manual’s.

Variable roles #

The time column’s scale is decided once, for the whole column, and a date-only value is anchored at local midnight. Choosing the scale per value would mix small ordinals with epoch milliseconds in one vector, and the spacing between a value on each scale is not a cadence – so one cell that parses only as a date makes the column a date column, and every numeric-looking cell in it is then a parse failure rather than a year or an index. Date.parse reads a bare YYYY-MM-DD as UTC midnight but a YYYY-MM-DD HH:MM as local time, so a column holding both would carry two anchors; the module appends T00:00:00 to a date-only string so that every value is local, which is the scale the plots’ time axis and the result tables render with. Fewer than three parseable rows leave nothing to space, so the run is refused with a toast.

Validation. On the airline table all 144 months parse as dates and no row is dropped; the Time axis rows read 144 observations, regular spacing, 0 gaps and 0 duplicates.

The seasonal period is the calendar cadence nearest the median spacing on a log scale, or the observations per day below about seventeen hours. The median spacing is robust to the few gaps and duplicates a real file carries, and the distance to a cadence is judged as a ratio because the cadences are themselves ratios apart – daily to weekly is 7×, weekly to monthly 4.3× – while the jitter inside one is small: months of 28 to 31 days differ from the 30.44-day cadence by under 8%, so they all read as monthly without a table of thresholds. At more than one observation a day the cycle is the day itself and the period is one over the spacing in days, rounded – 24 for hourly data, 3 for eight-hourly – with the cut at $\sqrt{0.5}$ of a day, about seventeen hours, the geometric midpoint between one and two observations a day. A numeric index carries no cadence, so it resolves to 1 unless overridden, and a yearly cadence resolves to 1 as well, since no shorter cycle is there to model.

Validation. The airline months are 28 to 31 days apart; the median snaps to monthly, the frequency line reads Monthly spacing – frequency = 12 (yearly cycle), and findfrequency confirms 12 from the series itself.

Spacing, gaps and duplicates are measured inside each group and the per-group medians pooled. The interleaved timeline of two groups sampled over the same span alternates between them, so its spacings are the offsets between the groups – near zero where their timestamps coincide, read as duplicates, and a jump at every transition, read as a gap – and not the cadence either group was sampled at. The gap and duplicate counts are therefore per group, the frequency is read off the median of the per-group median spacings, and the cross-group median is only the fallback for a run in which no group has two spaced observations.

One seasonal period is resolved for the whole run, and a warning names the groups’ own periods when they differ. The R runners take one frequency and build every series at it; fitting each group at its own period would make the seasonal terms, strengths and forecasts of one card incomparable across its groups and would need an override per group. So the resolved period is global, and on a date axis with more than one group each group’s median spacing is snapped to a period on its own – when the distinct results number more than one, every result card opens with a warning listing them beside the period used, and the remedy is to split the analysis by cadence.

Pre-flight (Check data) #

A gap is a spacing wider than 1.5 × the group’s median spacing, and a duplicate a spacing of zero. Calendar jitter is small – the longest month is 31/28 = 1.11 times the shortest – while one skipped observation doubles the spacing, so a factor of 1.5 separates the two with margin on either side; a spacing tighter than the median, two readings close together, is not a gap and is counted only when it is exactly zero. The JavaScript side counts for the dataset-wide rows and the R side for the per-group table, on the same rule.

Validation. The airline spacings run from 28 to 31 days against a median of 31, so the gap threshold is 46.5 days and no month is one: 0 gaps, 0 duplicates, regular spacing Yes.

Missing values are interpolated before the difference counts and the period search, and the raw series is kept for the outlier scan. ndiffs, nsdiffs and findfrequency take no NA, and stripping the blanks would shift every later value one step earlier and destroy the seasonal alignment the tests rely on, so a series with blanks is interpolated with forecast::na.interp – seasonal when the frequency allows – for those three computations only, and the raw series is used if the interpolation fails. ndiffs runs its default KPSS test at α = 0.05, up to two differences, on its own level and not the app’s significance setting; nsdiffs runs its seasonal-strength heuristic – one seasonal difference when the STL seasonal strength exceeds 0.64 – and only when the frequency is above 1 and the series holds at least two full cycles; both, and findfrequency, need eight non-missing values. findfrequency removes a linear trend, fits an autoregressive spectrum to what is left and reports the period at its peak when the peak is high enough to be one, else 1.

Validation. On the airline table the app suggests d = 1 and D = 1 – the series’ seasonal strength is 0.94 against the 0.64 cut – and detects the period 12.

Guerrero’s λ needs ten positive values and is estimated at the series’ seasonal period. forecast::BoxCox.lambda(method = "guerrero") is run on the series as a ts at the resolved frequency, so Guerrero’s coefficient-of-variation criterion compares the spread and level of whole seasons – the variation λ is meant to stabilise – and at frequency 1 those of consecutive blocks of two values; this is the reading forecast documents for a seasonal series and the one auto.arima(lambda = "auto") makes. Blanks are filled with forecast::na.interp for the estimate alone, since dropping them would slide every later value into the wrong season; a fill that leaves a zero or negative value reads as a failed estimate. Fewer than ten values leave the criterion nothing to fit, so the cell reads Not computed; a zero or negative value puts the series outside the power transform’s domain, and since the module does not shift a series the cell reads Not applicable; an estimate that errors reads Estimation failed. A series no longer than two seasonal periods reads Not estimated: at that length BoxCox.lambda returns λ = 1 without fitting anything, and printing that 1 as an estimate would report the identity as Guerrero’s answer. The same gates govern the transform stack’s automatic Box-Cox step, which skips such a series with the reason, and every training window the forecast’s holdout and rolling-origin validation refit the stack on – a window too short for the estimate is fitted untransformed, and the stack line counts those windows, since they score a different model than the one fitted to the whole series. A typed λ is applied whatever the length.

Validation. The airline table reads λ = −0.295 in the pre-flight – against 0.111 when the same values are read as a plain vector with no period – and 0.071 in a fit with Replace outliers on, since the stack estimates λ after the seven summer peaks were replaced.

Outliers are forecast::tsoutliers at the configured frequency, on the raw series, from ten non-missing values. On a seasonal series longer than two cycles a robust STL fit removes the seasonal component when its strength reaches 0.6, a super smoother fits the trend, and a residual below $Q_1-3\text{IQR}$ or above $Q_3+3\text{IQR}$ of the residuals is an outlier, in two passes with the flagged points interpolated between them; below ten values nothing is flagged. The scan runs on the raw series because it previews what the transform stack’s cleaner would replace, and the positions list is cut at twelve so a series of spikes does not turn the row into a wall of text.

Validation. On the airline table the seven points flagged are summer peaks – July of 1949, 1950 and 1951, August of 1958 and 1959, and July and August of 1960: the summer swing grows with the level, so an additive season leaves the largest peaks as residual extremes. They are genuine seasonal highs, not errors, which is the reader’s call to make before switching Replace outliers on.

The recommendation reads the seasonal difference first, then the regular one, and a short series is flagged beside either. A suggested seasonal difference implies a seasonal model whatever the regular count says; with none, a regular difference means a non-stationary series; a missing regular count means no test could run; and none of the three means a stationary series. The short flag – fewer than 30 non-missing values, about two cycles of monthly data with a buffer – is orthogonal to the recommendation on purpose: a short seasonal series is exactly the case where SARIMA and Holt-Winters absorb noise as signal, so the two lines appear together.

Transform stack #

Outlier replacement runs first, on the raw series at the configured frequency, and λ is then estimated on the cleaned series. forecast::tsclean is tsoutliers followed by interpolation of the flagged points and of every missing value; a point counts as replaced when its raw and cleaned values are both finite and differ by more than $\sqrt{\epsilon }$, the missing values it filled are counted apart, and a failure leaves the series untouched and is recorded rather than dropped. Running it before Box-Cox matters because a spike drags Guerrero’s λ: estimating λ on the cleaned series gives the power for the series the model will actually see.

Validation. On the airline table the stack replaces 7 points and then estimates λ = 0.071, where the pre-flight’s estimate on the raw series was −0.295.

λ is handed to the fitting function rather than applied beforehand, so the back-transform is mean-unbiased; the fits that take no λ transform locally. ets, auto.arima, Arima, naive, rwf, snaive and stlf accept lambda = with biasadj = TRUE, so their fitted values, forecasts and intervals return to the original scale corrected for the retransformation bias, while their likelihood-based numbers – information criteria, parameters, residuals – stay on the transformed scale, which is what the scale note on the card says. ma() and thetaf take no λ, and the periodogram and the breakpoint search need an already-transformed series, so those apply BoxCox by hand and invert it afterwards. A typed λ is one number for the run, used as typed for every series and group; blank estimates it per series. Below ten non-missing values, on a series with a zero or negative value – which the module does not shift – or when the estimate fails, the transform is skipped and the reason printed on the card, rather than an untransformed series being fitted in silence.

The stack carries no differencing step. ETS, theta and Holt-Winters model the level directly: a pre-differenced input hands them a series of increments whose level is noise, and their forecasts would have to be re-integrated by hand. ARIMA’s d performs the differencing inside its own likelihood, so it lives on that card, and the stack keeps only the shape transforms that compose cleanly with every method.

Analysis methods #

Exploration #

ACF and PACF run on the interpolated series at its configured frequency, through forecast::Acf and Pacf, with na.pass and a band at 1.96/√n. forecast::Acf extends the default lag range to two seasonal periods when that is longer than the usual 10·log₁₀(n), so the seasonal spike and its first echo stay on the correlogram – stats::acf would cut a monthly series off at lag 21; na.action = na.pass keeps a residual NA in place rather than closing the gap, since a stripped vector shifts every later lag. Lag 0 is omitted from the plot, as ggAcf omits it – it is 1 by construction and would only stretch the axis – and the band is the ±1.96/√n of the white-noise null, with n the observations used. Four non-missing values are the floor, below which the estimate is meaningless.

Validation. On the airline table the band is ±0.163 and the ACF runs to lag 24 – two periods rather than the 21 that 10·log₁₀(144) would give – with 0.948 at lag 1 and 0.760 at lag 12; the PACF’s −0.540 at lag 13 is the seasonal-difference signature.

The decomposition is STL with a periodic seasonal window and robust fitting, or stats::decompose; it needs a frequency above 1 and two full cycles, and the multiplicative form a positive series. s.window = "periodic" fixes the seasonal shape across cycles, which is what makes the seasonal strength read the configured season and nothing else, and robust = TRUE down-weights an outlying month so a spike lands in the remainder rather than in the trend; the classical forms are R’s own, the trend a centred moving average of one period. Fewer than two cycles leave no seasonal component to estimate – stl refuses outright – and a zero or negative value has no multiplicative share, so those runs are skipped with a coded note rather than failing inside R. A series with missing values is decomposed on the na.interp fill, which itself uses STL when the season allows and avoids the spectral artefacts of a plain linear fill; when that fill fails, a linear interpolation with end-value extension is the fallback, since stl takes no NA.

An anomaly point is a remainder more than three MAD from the remainder’s median, at an observed position. The remainder of a good decomposition is noise around zero, so its median and MAD – the constant 1.4826 scaling the MAD to a standard deviation – are a robust centre and spread that a few excursions cannot inflate, and three of them is the same cut the forecasting card uses on residuals; at least four finite remainders are needed for the spread to mean anything. A position whose value was interpolated carries the fill’s own fit, not data, so an excursion there is excluded as an artefact of the interpolation.

Validation. On the airline table the STL remainder’s excursions are the later summers, whose swing outgrows the fixed periodic season: 18 points are flagged, June to August of 1956 to 1960 among them.

Trend and seasonal strength are Hyndman’s $1-\mathrm{Var}(R)/\mathrm{Var}(T+R)$ and $1-\mathrm{Var}(R)/\mathrm{Var}(S+R)$, floored at 0, on logged components for a multiplicative decomposition. The definitions are those of the tsfeatures package’s stl_features, so the numbers are comparable with the literature that uses them; the floor at 0 catches a remainder larger than the component-plus-remainder sum, which the formula would otherwise print as a negative strength. A multiplicative decomposition is additive in the log, and its remainder oscillates around 1 whatever the noise level, so the additive formula applied to the raw multiplicative parts returns a number pinned near 1 regardless of the data – logging the three components first restores the meaning.

Validation. On the airline table the STL reads trend 0.942 and seasonal 0.618; the pre-flight’s nsdiffs reads a seasonal strength of 0.94 on the same series because its own STL uses a different seasonal window, which is why the two numbers are not expected to agree.

The stationarity tests are tseries::adf.test, kpss.test with the selected null and pp.test in its Z(t_alpha) form, on the interpolated series, from eight non-missing values; a p outside the table prints as a bound. adf.test fits drift and trend and chooses its lag order as ⌊(n − 1)^(1/3)⌋, pp.test the same alternative with a Newey-West correction and its Bartlett truncation lag from the series length, and the Z(t_alpha) form is the t-ratio statistic, the only PP variant on ADF’s scale – hence PP (Zτ) in the table; kpss.test takes null = "Level" or "Trend" from the selector. All three read their p-values off interpolation tables and return the table’s endpoint when the statistic falls outside it, announcing the truncation only in a warning, so the runner captures the warning and the cell prints < 0.01 or > 0.10 rather than the endpoint as if it were exact; a bounded cell keeps its significance class and stops re-rendering as a number. Unit-root tests take no na.action, so the series must be complete: the interpolated copy is used, and when the interpolation itself fails the tests are skipped rather than run on a stripped vector, whose shifted index would bias every statistic. Below eight values adf.test’s lag order is 1 on a handful of differences, so the floor is shared by every test.

Validation. On the airline table ADF is −7.32 with lag 5 and PP −5.05 with lag 4, both < 0.01, and the level-form KPSS 2.74 with lag 4, also < 0.01 – every cell a bound, since the statistics lie far beyond the tables.

KPSS is judged at the assumption level and the unit-root tests at the main level, and the conclusion flips the null’s own claim. KPSS holds “the series is stationary” – the data-are-fine null every assumption check holds – so it reads settings.assumptionAlpha like the normality and sphericity tests; ADF and PP hold the opposite null, so they keep the main significance level, which is what keeps the pair moving the same way when the user parts the two levels. A rejection on a stationarity null is a non-stationary verdict and a rejection on a unit-root null a stationary one, so the conclusion column is computed from the null code and the row’s own α, and the α printed in each conclusion is the one that row was judged at.

The joint reading pairs the unit-root verdict with KPSS in a 2 × 2 and refuses to pair when ADF and PP disagree. A single test can only fail to reject, so only the two opposite nulls together tell a stationary series from a sample too weak to decide: reject-and-not is stationary – trend-stationary when the KPSS null was trend, which points at detrending rather than differencing – not-and-reject a unit root, both a series neither description fits, neither an uninformative sample. When both unit-root tests ran and one rejects while the other does not, there is no unit-root verdict to pair, so the line says so rather than picking ADF over PP.

Validation. On the airline table ADF and PP reject and the level-form KPSS rejects too, so the line reads the both-reject case – a deterministic trend, a break or changing variance – which is the series’ growing swing around its trend.

The raw Ljung-Box uses Box.test at max(10, 2m) lags, capped at n/5. Hyndman’s rule for a seasonal series is 2m lags so the seasonal echo is inside the window, and 10 otherwise; the cap at a fifth of the series keeps the statistic’s χ² approximation sound on a short series. Its null is no autocorrelation, so it reads the assumption level.

Validation. On the airline table the lag is min(24, 28) = 24 and Q(24) = 1606, p < .001.

The seasonal-period search is a stepwise Fisher’s g-test on an untapered, unpadded periodogram, with a Bonferroni bound for the p, a ±2-bin flattening per step, at most ten steps and no second adjustment. spec.pgram runs with taper = 0 and fast = FALSE because both break the i.i.d. exponential null Fisher’s g assumes – a taper pulls the ordinate degrees of freedom below 2, and zero-padding inflates the ordinate count while correlating neighbours; demean and detrend are on so the zero-frequency mass does not pass for a cycle. The p is the bound $N{(1-g)}^{N-1}$ on Fisher’s exact alternating sum, which is numerically unstable for the long series the module may see (N ≈ 2500) and errs conservatively. At each step the largest remaining ordinate is tested, then its main lobe – two bins either side – is flattened to the periodogram’s median, since a cycle that is not an exact Fourier frequency smears across its neighbours and a single-bin removal leaves shoulders the next step rediscovers as a slightly shifted peak; the loop stops at ten steps or when a step’s p reaches 0.5, since the residual periodogram is then white-noise-like. Survivors are those with p at or below the significance level and a period between 2 observations and half the record; the g of a later row is measured on the flattened working periodogram, so it exceeds the row’s share of the untouched total, which is what the Share of variance column prints. The bound already controls the family-wise error over every ordinate, so applying the global p-value adjustment over the K peaks would correct the same family twice.

Validation. On the airline table 10 peaks are tested and 7 survive at α = 0.05: 12 months with g = 0.502 and p = 2 × 10⁻²⁰, then 6, 4, 3 and 2.4 months, and two weak long periods, 8 and 9.6, at p = 0.029 and 0.044. The 6-month row’s g is 0.404 on the flattened periodogram against a 15.3% share of the untouched one.

A row is a harmonic when its Fourier-bin interval overlaps a stronger, longer row’s interval divided by an integer of two or more. A non-sinusoidal cycle of period P puts real peaks at P/2, P/3, …, which the test rightly calls significant; the spectral resolution already reported as the period’s bounds is the tolerance the ratio deserves, so the check is an interval overlap rather than a rounding of P/p, and the strongest matching parent wins since the candidates are walked in g order. The row keeps its place in the table – it is a real peak – and is muted and annotated rather than dropped.

Validation. On the airline table 6, 4, 3 and 2.4 months are marked ×2, ×3, ×4 and ×5 of 12; 8 months is not, since 12/8 rounds to 2 and the 6-month interval [5.9, 6.1] does not reach the 8-month one [7.8, 8.2].

Smoothing #

A moving average is forecast::ma(order = k, centre = TRUE), seeded from the seasonal period, floored at 2, and under Box-Cox transformed by hand and back on the median scale. ma applies the 2×m weighted form at an even order, so a window of 12 is the classical centred seasonal trend estimator and the curve spans m + 1 points; the seed is the active frequency because the window that removes a season is one season wide, and it falls back to 7 – a week of daily data – when no period is available. A window of 1 is the series itself and 0 is undefined, hence the floor. ma accepts no lambda, so the series is BoxCox-transformed before the filter and InvBoxCox-inverted after; a filter has no residual variance to feed the bias-adjustment formula, so the inverted curve is a median rather than a mean and the card says so. Fewer non-missing values than the window leave no complete window, so the fit is skipped.

Validation. On the airline table a window of 12 – the seed at frequency 12 – smooths 132 of 144 points, six lost at each end.

SES, Holt and Holt-Winters are forecast::ses, holt and hw with biasadj = TRUE, on an interpolated series, from 4, 5 and 2m non-missing values. The three are the ETS models ANN, AAN or AAdN, and AAA / MAM with their damped forms, and ets needs a handful of observations more than parameters – four for a level, five for a level and a trend, two full cycles for a season – below which it refuses to fit or returns a degenerate one. They cannot see through a missing value – they either refuse or collapse every smoothing parameter to its lower bound – so a series with blanks is filled with na.interp before the fit and the fill is counted on the transform-stack line; a moving average reports its NA windows honestly and is left alone. hw rejects a seasonal period above 24 outright and multiplicative seasonality on a Box-Cox-transformed series, and multiplicative indices are undefined on a series touching zero, so each of those is a coded skip before R is asked. The sigma reported is $\sqrt{{\sigma }^2}$ of the fitted model and the parameter table is model$par split into the four smoothing weights and everything after them, which ets stores as the initial states.

Validation. On the airline table SES estimates α = 0.9999 – a naïve forecast, RMSE 33.6 – and multiplicative Holt-Winters ETS(M,A,M) with α = 0.315, β = 0.007, γ = 0.598, σ = 0.041, RMSE 10.6 and MAPE 2.91%.

The residual diagnostics are checkresiduals’ Ljung-Box – 2m or 10 lags, capped at n/5, with fitdf the model’s degrees of freedom – a residual ACF banded at the confidence level, and a Shapiro-Wilk on at most 5,000 seeded residuals. Box.test subtracts fitdf from the lag to get its degrees of freedom and reports df ≤ 0 rather than erroring, so a lag that cannot clear forecast::modeldf – the number of parameters the fit estimated, the constant included on an ARIMA – skips the test with a note instead of printing a meaningless df. Residual NAs reach the tests un-stripped, since closing the gaps shifts every lag estimate. shapiro.test caps at 5,000 observations, so a longer residual vector is judged on a draw of 5,000 under settings.randomSeed, and the draw is reported beside the verdict; three residuals are its floor. Both tests hold the fitted model as their null – white noise, normal residuals – so they read the assumption level.

Validation. On the airline table the multiplicative Holt-Winters residuals give Q(24) = 40.2 on 24 − 16 = 8 degrees of freedom, p = 3 × 10⁻⁶, and W = 0.990, p = 0.43; the automatic ARIMA’s give Q(24) = 41.8 on 23 degrees of freedom, p = 0.010, and W = 0.978, p = 0.021.

ARIMA and SARIMA modelling #

Auto mode is forecast::auto.arima under the chosen criterion with stepwise and approximation both on by default and both off under Exhaustive search; ten non-missing values are the floor. The stepwise search with an approximated criterion is Hyndman and Khandakar’s default and finishes in seconds on a long series; scoring every order exactly regularly finds a lower criterion and costs minutes on a long or high-frequency series, so it is a switch rather than the default. The order is read back through arimaorder, which reports P, D and Q only on a seasonal fit, and the constant through the coefficient names – drift on a differenced series, intercept on an undifferenced one – so the heading can say with drift the way summary.Arima does. Below ten values auto.arima has too few differences to choose an order on, and a manual fit’s neighbouring grid is meaningless.

Validation. On the airline table the default search returns ARIMA(0,1,1)(0,1,0)[12] with ma1 = −0.318, AICc 1020.7, RMSE 11.2 and MASE 0.258.

Manual mode is forecast::Arima with include.constant, its seasonal orders forced to zero at frequency 1, and a 3-wide neighbouring grid per axis with d and D held fixed. Arima estimates a seasonal term at period 1 without complaint – confounded with the non-seasonal one and reported back by arimaorder as if it were not there – so the seasonal row is disabled in the panel and zeroed in the runner when the frequency is 1. include.constant fits a mean on an undifferenced series and a drift on a once-differenced one, as auto.arima does, and nothing on higher differencing. The grid varies p, q, P and Q one step either way of the chosen order – {0, 1, 2} at a boundary of 0, so the comparison stays symmetric and shows whether moving up to 2 helps – and holds d and D fixed, since information criteria across different differencing are not comparable; that is 9 candidates on a non-seasonal frequency and 81 on a seasonal one, each fitted in full, which is why the grid runs only in manual mode and carries its own progress tick. The bounds – p, q ≤ 10; d ≤ 3; P, Q ≤ 5; D ≤ 2 – are the widget’s min/max, and a blank or out-of-range order stops the run with a toast before anything is fitted – the option builder’s fallback to 0 serves only the pristine panel.

Validation. On the airline table the airline model ARIMA(0,1,1)(0,1,1)[12] fits ma1 = −0.309 and sma1 = −0.107 at AICc 1021.2; 79 of its 81 neighbours fit, the best of them ARIMA(1,1,2)(1,1,2)[12] at AICc 1013.2.

Coefficient intervals are estimate ± z × SE at the global confidence level, with the p from the normal; the training-set accuracy is forecast::accuracy on the back-transformed fitted values; and above period 24 the card warns. Arima reports var.coef, whose diagonal gives the standard errors, and the asymptotic distribution of the estimates is normal, so a z rather than a t is the right reference. accuracy() scores fitted(), which biasadj = TRUE returns on the original scale, so under Box-Cox the four accuracy rows are the only original-scale numbers on the card while the criteria, σ² and residuals stay transformed – the scale note. The seasonal period past which auto.arima quietly drops to seasonal differencing alone and a manual seasonal term stops converging is about 24 – Hyndman’s own ceiling for SARIMA – so the card warns there and names dynamic harmonic regression and STL-based methods as the alternative.

Forecasting #

The seven forecasters are naive, rwf(drift = TRUE), snaive, forecast(ets(…)), forecast(auto.arima(…)), a theta wrapper and stlf(robust = TRUE, s.window = "periodic"), each with lambda and biasadj = TRUE where it takes them, and each with a coded precondition. The set covers the baselines, the two automatic model families, the M3-winning theta method and the STL route for long or awkward seasons, with the prediction intervals asked at 80 and 95. stl needs two whole periods, snaive one and both a frequency above 1, and thetaf’s ACF step fails on any missing value the fill could not replace; each of those reports only an opaque error from R, so the runner checks the precondition first and the row reads Not applicable with the reason, while a fit that breaks on the data reads Failed. thetaf takes no lambda, so its input is Box-Cox-transformed by hand and the mean and fitted values are inverted with biasadj and a forecast variance recovered from the 95% interval’s half-width, so the theta forecast is mean-unbiased like the others on the same plot; and thetaf decomposes multiplicatively whenever its own seasonality test fires, which is undefined on a series crossing zero, so that test is mirrored and a non-positive seasonal series is adjusted additively instead.

The holdout is the last min(h, ⌊0.2 n⌋) observations of a series with 30 non-missing values and a training slice of 10, scored on its finite actuals, with the transform stack refitted on the training slice. A fifth of the series is the usual ceiling on a validation window, and the horizon is what the forecast will actually be asked for, so the smaller of the two is scored; below 30 observations any window worth scoring leaves too little to fit on, so the block is skipped and the validation line names the gate. The window is positional while the gate counts non-missing values, so the errors are taken on the finite actuals inside it. The stack is refitted on the training slice because an outlier replacement or a λ estimated with the window in view would leak the window into the fit it is meant to test.

Validation. On the airline table the holdout is min(12, 28) = 12 months, the year 1960; Auto-ARIMA scores RMSE 23.9 and MASE 0.61 on it.

Rolling-origin cross-validation places its origins a full horizon apart, runs min(requested, ⌊(n − 10)/h⌋) of them, and scores only the origins every fitted method reached. Disjoint windows keep any observation from being scored twice, so the pooled error is an honest average rather than an over-weighted recent stretch; each origin therefore costs a horizon of length, and the count the series can support is what the runner uses, reported as a of n when it is fewer than asked. A method that fails at some origin would otherwise leave that origin in the pool for its rivals and out for itself, so an origin scores only when every method in the race forecast it, and a method that reached no origin at all leaves the cross-validation rather than emptying the shared set. The maximum of 20 origins bounds the cost – 20 × 7 fits per series – on the headless path as well as in the widget.

Validation. On the airline table five origins a year apart all score: Auto-ARIMA RMSE 27.0 and MASE 0.71, ETS 30.5 and 0.81, Theta 33.3 and 0.90, seasonal naïve 42.5 and 1.33, STL + ETS 60.8, Drift 80.8, Naive 90.4 and 2.30.

MASE divides by the in-sample seasonal-naïve MAE of the slice the fit trained on, MAPE is undefined unless every scored actual is positive, and one metric ranks the whole series – CV RMSE, else holdout RMSE. Hyndman and Koehler’s scaling uses the training-set naïve error at lag m – lag 1 on a non-seasonal frequency – so a MASE below 1 means better than repeating last season on data the method also saw; a zero denominator leaves it undefined. MAPE divides by the actual, so an actual at or near zero – within 10⁻⁶ of the largest actual – makes it undefined rather than enormous. CV and holdout score different data, so ranking on one metric for the whole series is the only comparison that holds: CV RMSE whenever at least one origin scored, holdout RMSE otherwise, a method missing the metric sorting after every scored one and a method that never fitted sorting last whatever the metric holds; the rank header names the metric so the order is never read as a cross-validated one when it is not.

In-sample anomalies are residuals more than three MAD from their median, against the raw series, from five finite residuals. The residuals score the raw observations rather than the transformed series, so a point tsclean replaced stays eligible for a flag – the flag says the model saw the observation as unusual, whatever the stack did to it – and the 3·MAD rule is the exploration card’s, so the two markers agree on what an excursion is.

Future timestamps step in calendar months when the spacing is a whole number of months within tolerance, from the anchor rather than cumulatively; every other cadence extends linearly. A monthly spacing is 28 to 31 days and a constant-millisecond extrapolation drifts off the month boundaries within a year, so a spacing that is a whole number of 30.44-day months – within 3% of it, floored at two days, which absorbs month-length jitter at one month while still rejecting a 28-day four-weekly cadence – is walked in calendar months. Each step is measured from the anchor and clamped to the target month’s length, because a cumulative walk overflows on an end-of-month anchor – 31 January plus one month is 3 March – and skips a period outright. The prediction intervals are fixed at 80% and 95%, the two levels every forecast method returns by default and the two the plot and table can show without a level control of their own.

Spectral / change-points #

The periodogram is spec.pgram with a 10% taper, detrended, on the true Fourier grid, smoothed for the plot by a modified Daniell kernel of span from √n, with the five strongest local maxima up to period n/2 tested by Fisher’s g on a matching untapered periodogram. fast = FALSE keeps the true Fourier grid: zero-padding to a highly composite length invents bins below 1/n, periods longer than the record itself. The taper is the standard 10% split cosine, which reduces leakage from strong peaks into their neighbours; the raw periodogram carries about two degrees of freedom per ordinate and reads as a noise comb, so a modified Daniell smooth of span max(3, 2⌊√n/4⌋ + 1) is drawn beside it – beside, not instead, since Fisher’s g is defined on the raw ordinates – and its bandwidth is stated under the plot. Peaks are local maxima, so adjacent bins of one peak do not fill the list; one or two observed cycles are not evidence of a period, so the ranked list stops at n/2 while the plot keeps the whole spectrum. The g of each peak is taken on a second, untapered periodogram, because a taper pulls the ordinate degrees of freedom below 2 and breaks the i.i.d. exponential null; here it is the peak’s share of the whole spectrum with the same Bonferroni-bound p as the exploration card, so the first rows of the two cards agree and later rows differ, the exploration card having flattened each stronger peak before testing the next. Detrending a constant series leaves an identically zero spectrum in which every bin ties, so that case is skipped rather than ranked; eight observed values are the floor for a spectrum with any resolution.

Validation. On the airline table the 72 ordinates are smoothed at span 7, bandwidth 0.0125 cycles per observation; the peaks are 12 months (g = 0.502, p = 2 × 10⁻²⁰), 6 (g = 0.153, p = 0.0005), 4, 3 and 48 (p = 1), where the exploration card’s 6-month row reads g = 0.404 on its flattened periodogram.

Breakpoints are strucchange::breakpoints(y ~ 1) or (y ~ t) with h the minimum segment fraction, the number of breaks by BIC, and confint intervals clamped to the record; segment means are taken on the untransformed series, and twenty observed values are the floor. y ~ 1 is the pure mean-shift model and reads any trend as a staircase of level steps; y ~ t is the segmented-trend form. breakpoints selects the number of breaks by BIC over the dynamic-programming search, and h caps how many can be placed: the widget’s 5 to 45% keeps room for at least two segments, since the search caps out well before half the record, and the resulting length in observations is printed so a clamped value is visible. confint errors outright when no break was found, hence the separate no-break branch; its interval can run past either end of the record and is clamped there. A mean taken on the Box-Cox scale and pushed back through InvBoxCox is a median, not the mean of the data the plot draws, so the segment means are read off the untransformed interpolated series. Below twenty values the segments at any allowed h hold too few observations to estimate a level or a line. strucchange drags zoo, sandwich and lattice in with it, so it is loaded only when the search is on, and a load failure degrades to a warning on the breakpoint section while the periodogram still runs.

Validation. On the airline table the mean-shift model at h = 15% – a minimum segment of 21 observations – places five breaks after observations 26, 50, 77, 98 and 123, the first dated within [24, 28], with segment means 134, 188, 239, 318, 379 and 464; the level-and-trend model places none.

The length gates count observed values, and the interpolation is counted on the transform-stack line. na.interp fills every gap once two real values exist, so a 60-row series holding four real values would pass an eight-value gate on its filled length; the gates read the pre-imputation count instead, the fills join the transform-stack record rather than passing for observations, and the closing line names both counts.

Imputation, missing data, and dropped rows #

A blank series value is kept in place and filled per card, never by dropping its row, and the global missing-data pass is left to act upstream. ts() reads no time column – position is time – so deleting an incomplete row moves every later observation one period earlier and misaligns the seasons, while a kept NA keeps the grid and lets each fit decide: Arima and auto.arima take a gap through their Kalman filter, the lag-walk forecasters naive, rwf and snaive step over it, stlf interpolates inside mstl on its own, and the fits and tests that take no NA – ses, holt, hw, the unit-root and KPSS tests, Box.test, stl, spec.pgram, breakpoints – are handed an na.interp fill, which follows the series’ trend and season, and the fill is counted on the card. The forecasting card’s ets and thetaf take no NA either and get the same fill, while the forecasters beside them keep fitting through the gap; the fill is redone on every training slice – the full series, the holdout’s and each cross-validation origin’s – so no value past a slice’s end shapes it, and accuracy is scored against the observed values alone, so a filled point never scores. The only rows the module drops are those with no place on the axis – an unreadable time, or a blank group under a grouping variable – and their count is reported on the pre-flight and in every card’s warning. The global missing-data pass runs on the dataset before any module reads it, so the module cannot switch it off and does not claim to: listwise deletion removes the incomplete rows before the time axis is read, and they reappear as gaps in the spacing diagnostics rather than as missing values; single-value imputation fills every blank in a column with one constant, the flattest fill a seasonal series can take, and leaves nothing for a card to disclose. The manual describes both and recommends pairwise deletion, under which the blanks reach the cards as NA and each card’s own fill is the one applied.

Validation. On the airline table with three months blanked – May 1953, November 1956, February 1959 – the forecasting card races all seven methods, ETS and Theta on the filled series (3 values filled), and ranks Auto-ARIMA first at CV RMSE 27.4, ahead of ETS at 30.4 and Theta at 33.3; under listwise deletion the pre-flight reads 141 observations, 3 gaps and 0 missing values.

Reproducibility #

Citations follow the run, not the repertoire: forecast and its paper always, tseries only when a stationarity test was selected, strucchange only when the break search ran, and each method from the branch that computed it. A reference list naming every method the module offers would be useless as a Method paragraph, so each runner cites what it did. The pre-flight cites KPSS because ndiffs decides by one, Box-Cox and Guerrero because it estimates λ, the super smoother because tsoutliers fits one, and STL and the strength measure only at a frequency above 1. Exploration cites the ACF and PACF, the Ljung-Box test, the periodogram, Fisher’s g and the Bonferroni bound behind its seasonal-period search whatever the options, each stationarity test only when it was ticked, the decomposition method and the 3·MAD rule only when a decomposition ran. The fitted cards cite their model family – the smoother chosen, ARIMA and the automatic search in auto mode, the seven forecasters – with the three information criteria and MASE where the card prints them, each diagnostic only when it was on, rolling-origin cross-validation and the anomaly rule only when they ran; and the transform stack cites the super smoother and STL only with Replace outliers on, Box-Cox only with the override on and Guerrero only when λ was estimated rather than typed. The one seeded draw – the Shapiro-Wilk sub-sample above 5,000 residuals – is argued with the residual diagnostics; nothing else in the module draws a random number.