Time Series in R Exercises: 50 Real-World Practice Problems

Explanation: frequency is the number of observations in one natural cycle, so 4 means quarterly and 12 means monthly. The start argument takes c(year, period_within_year), not a date. Once R knows the frequency it prints the quarterly grid for free and every downstream function (window, aggregate, decompose) inherits that calendar. A common mistake is passing start = 2019 alone, which silently assumes period 1 and works here but breaks the moment your series starts mid-year.

Explanation: The blanks in 2021 and 2022 are not missing values; the print method simply lays every series out on a full January-to-December grid, so partial first and last years look ragged. This is the single most common source of off-by-one bugs in R time series work: if you had written start = c(2021, 1), every observation would sit six months early and your seasonal model would learn the wrong month effects. Check start() and end() after constructing anything.

Explanation: start() and end() return c(year, period) pairs, which is why the [1] and [2] indexing is needed. The arithmetic should tie out: 12 years times 12 months is 144 observations, so there are no gaps. If length() disagrees with what the start and end dates imply, you have a ragged or irregular series and every seasonal method will misbehave. tsp(AirPassengers) gives the same information in one shot as c(start_time, end_time, frequency) in decimal years.

Explanation: AirPassengers[73:108] returns the same 36 numbers but as a plain vector with no frequency, so decompose() and forecast() would then refuse to run or would treat the data as annual. window() preserves the tsp attribute, which is why the result still prints as a month grid. This is also the idiomatic way to carve a training set out of a series: window(x, end = c(1958, 12)) for train, window(x, start = c(1959, 1)) for test.

Explanation: time() returns decimal years (1959.000, 1959.083, ...) so flooring it yields the calendar year, while cycle() returns the position within the cycle, which for a monthly series is the month number. as.numeric() on the series strips the ts attributes so the tibble stores plain doubles. Going ts to tibble like this is the standard bridge into dplyr and ggplot2 work; going back the other way is just ts(df$ppm, start = c(1959, 1), frequency = 12).

Explanation: An mts object is stored as a matrix with a tsp attribute, which is why dim() and colnames() work alongside frequency(). The frequency of 260 encodes trading days per year rather than calendar days, a convention worth checking before you trust any seasonal method on financial data. To pull one index out and keep it a series, use EuStockMarkets[, "DAX"], and to bind separate series back into an mts, use ts.union(), which aligns on the calendar and pads with NA instead of recycling.

Explanation: A ts object has no dates, only a start and a fixed step, so it assumes every consecutive row is one period apart. Feeding these 7 rows straight into ts(readings$kwh, frequency = 7) would silently compress 10 calendar days into 7 slots and shift every later observation. Densifying against a complete date spine first is the only safe conversion path for irregular data. Once the gaps are explicit you can decide what to do with them: zoo::na.approx() to interpolate, tidyr::fill() to carry the last value forward, or leave the NA and use a model that tolerates it.

Explanation: Duplicated timestamps and dropped intervals are the two failure modes that quietly corrupt every downstream lag, difference and seasonal index, because both shift observations into the wrong period. Doing the audit as a one-row summary means you can bind it across many feeds and monitor the counts over time. In production you would usually promote this to a hard assertion, for example stopifnot(n_distinct(events$ts_min) == nrow(events)), so a bad load fails loudly instead of producing a plausible but wrong forecast.

Explanation: aggregate.ts() defaults to nfrequency = 1, which is why summing gives calendar-year totals and the result prints with Frequency = 1. Aggregation is not neutral: the annual series is smooth and trending because summing has destroyed the seasonal pattern entirely, which is exactly what you want for a capacity plan and exactly what you do not want if the question is about peak months. Note this only works cleanly because the series contains whole years; a partial final year would be summed as though it were complete.

Explanation: nfrequency must divide the original frequency evenly, so 12 to 4 and 12 to 1 are legal but 12 to 5 is not. Use FUN = mean for level-type quantities such as concentrations, prices or temperatures, and FUN = sum for flow-type quantities such as units sold or passengers carried. Mixing those up produces numbers that are wrong by a factor of three at quarterly resolution and are easy to miss because they still look plausible.

Explanation: Differencing costs you one observation, which is why the differenced series starts in 1872 while Nile starts in 1871. That bookkeeping matters when you line a differenced series up against a predictor: lengths no longer match and R will recycle rather than error. Differencing is the standard fix for a wandering mean, and the d in ARIMA(p, d, q) is literally the number of times diff() was applied before the model was fitted.

Explanation: The difference of logs approximates a proportional change, so 0.0522 reads as roughly 5.2 percent growth from January to February 1949. That approximation stays close for changes under about 10 percent and drifts beyond it, where exp(d) - 1 is the exact figure. Logging first is what turns a multiplicative seasonal pattern into an additive one, which is why the airline model in Section 5 is fitted to log(AirPassengers) rather than the raw counts. Logs require strictly positive data; for series containing zeros, log1p() or a Box-Cox transform is the usual substitute.

Explanation: dplyr::lag() shifts values down one row and pads the top with NA, which is why January has no growth rate. Two traps are worth knowing. First, stats::lag() does something completely different: it shifts the ts index rather than the values, so loading dplyr after stats is what makes this code behave as expected. Second, lag() respects row order but not calendar gaps, so always arrange() by date first and confirm there are no missing periods, otherwise a skipped month silently becomes a two-month change.

Explanation: sides = 2 centres the window on each observation, so five months at each end of the series have no full window and come back NA. That edge loss is unavoidable and is the reason a trailing moving average (sides = 1) is preferred for live dashboards: it costs no data at the recent end, at the price of lagging turning points by half the window. Note the qualifier stats:: is needed because dplyr also exports a filter(), and calling the wrong one here throws a confusing error about an unused argument.

Explanation: Because the lag equals the frequency, every value here compares a month to the same month in the prior year, so the entire seasonal shape cancels and what remains is year-over-year growth. Seasonal differencing costs a full cycle, hence the result starting in 1950. Note diff(x, lag = 12) differences once at lag 12, whereas diff(x, differences = 12) differences twelve times at lag 1; the argument names look similar and mixing them up destroys the series.

Explanation: Rebasing makes series with different units comparable on one axis, which is why price indices, stock charts and benchmark comparisons nearly always use it. The as.numeric() matters: Nile[1] already returns a bare number here, but on a multivariate series the equivalent subsetting can drag ts attributes along and produce a length mismatch. Watch the base-period choice, because an unusually high or low first observation makes every subsequent value look permanently depressed or inflated; a common alternative is to divide by the mean of the first full cycle.

Explanation: plot() dispatches to plot.ts() when handed a series, which is why the x axis comes out in calendar years without any effort. Read this chart before you fit anything: it shows an upward trend, a strong annual cycle, and crucially a seasonal swing that widens as the level rises. That widening is the visual signature of multiplicative seasonality, and it is the reason later exercises log the series first. Getting into the habit of plotting before modelling catches structural breaks and outliers that no automatic order-selection routine will flag for you.

Explanation: A standard time plot mixes trend and season together, so it is hard to tell whether the seasonal shape itself is changing. monthplot() separates them: the slope within each panel is the trend for that month, and the difference between the panel means is the seasonal effect. Here every month climbs at a similar rate, which says the seasonal pattern is stable and additive rather than growing, so decompose(co2, type = "additive") is the right choice. Despite the name, monthplot() works for any frequency, including quarterly and weekly series.

Explanation: cycle() turns the ts index into a month factor, which is what makes the formula interface work here. The boxes do double duty: their centres trace the seasonal profile, and their heights reveal that the peak months vary far more in absolute terms than the quiet ones. That heteroscedasticity is a second argument for a log transform, since logging compresses the large values and makes the spread roughly equal across months. Be aware the trend inflates every box, so this view exaggerates spread compared with boxplots of a detrended series.

Explanation: In a multiplicative decomposition the factors are ratios centred on 1, so 1.2266 for July means July traffic runs about 23 percent above the trend while 0.8012 for November runs about 20 percent below it. The $seasonal element repeats these twelve numbers across all 144 months; $figure is the compact one-cycle version. The factors are constant by construction, which is the main limitation of classical decomposition: if the seasonal shape genuinely evolves, use stl() with a finite s.window instead.

Explanation: STL fits trend and seasonality with repeated loess smoothing, which gives it two advantages over decompose(): it is robust to outliers, and it produces trend estimates right to the ends of the series instead of NA padding. s.window = "periodic" forces one fixed seasonal shape, equivalent in spirit to classical decomposition; setting it to an odd number such as 13 lets the shape drift slowly across the years. STL is additive only, so for a multiplicative series you decompose log(x) and exponentiate the components afterwards.

Explanation: Compare the adjusted July 1949 figure of 120.7 with the raw 148: the summer spike was mostly the calendar, not underlying growth. Getting the arithmetic right is the whole exercise, since dividing an additive decomposition or subtracting a multiplicative one leaves a residual seasonal pattern that is easy to miss. Note the adjusted series is not smooth; it still carries the remainder, which is exactly right because seasonal adjustment removes the calendar effect and nothing else. Official agencies use X-13ARIMA-SEATS rather than this simple ratio, mainly to handle trading-day and holiday effects.

Explanation: The amplitude climbs from 44 to 232 while the level roughly quadruples, so the swing grows in step with the level rather than staying constant. That is precisely the definition of multiplicative seasonality, and it rules out an additive decomposition, which assumes one fixed seasonal offset for every year. The practical consequence is that you either fit a multiplicative model or you log the series and fit an additive one, which are close to equivalent. As a quick check, the amplitude as a fraction of the yearly peak stays near 0.3 to 0.4 throughout, whereas a truly additive series would show that ratio shrinking.

Explanation: These are the strength measures from Hyndman and Athanasopoulos, and both live on a 0 to 1 scale where values near 1 mean the component dominates the noise. The max(0, ...) guard matters because on a series with no real trend the ratio can exceed 1 and produce a meaningless negative score. Reading the numbers here: co2 is almost pure trend with a strong stable season, so a seasonal model is mandatory and a plain random walk would be badly wrong. At scale these features are what let an automated pipeline route thousands of series to the right model family without a human inspecting each plot.

Explanation: Lag 0 is always exactly 1 because a series correlates perfectly with itself. The pattern here decays gradually rather than cutting off sharply, which points at an autoregressive component rather than a moving average one. With 100 observations the approximate significance threshold is 2 divided by the square root of 100, so 0.20; on that basis lags 1 through 8 are meaningful and lags 9 and 10 are noise. Extract the raw numbers with ex_4_1$acf when you need them programmatically, and remember the ACF assumes stationarity, so a slowly decaying ACF often just means the series needs differencing first.

Explanation: The PACF measures the correlation at each lag after removing everything the shorter lags already explain, so an AR(p) process shows large partial autocorrelations up to lag p and near-zero values after. Here lag 1 at 0.960 and lag 2 at -0.267 stand out against the 0.20 threshold while the rest fade, suggesting an AR(2) on the raw scale. The near-1.0 first partial is also a classic unit-root signature, which is why auto.arima() ends up differencing this series in Exercise 5.6. PACF starts at lag 1, unlike the ACF, because the lag-0 partial is undefined.

Explanation: The null hypothesis is that the first 10 autocorrelations are jointly zero, meaning the series is white noise. A p-value of 1.255e-14 rejects that decisively, confirming there is real structure worth modelling. Always pick type = "Ljung-Box" over the default Box-Pierce, since it has better small-sample behaviour. The same test is the workhorse for residual diagnostics later, but there you must pass fitdf equal to the number of estimated ARMA coefficients so the degrees of freedom account for the fitting, as Exercise 6.7 does.

Explanation: The direction of this test trips people up constantly. The null is that a unit root is present, so the series is non-stationary, and the alternative is stationarity. A p-value of 0.3107 therefore fails to reject non-stationarity, and the series needs differencing. Failing to reject is not proof of a unit root, only absence of evidence against it, which is why practitioners pair ADF with a KPSS test that flips the hypotheses. Exercise 4.6 shows that a single difference is still not enough here, a result you would never guess from the level plot alone.

Explanation: ndiffs() repeatedly applies a unit-root test until it stops rejecting, and nsdiffs() does the seasonal equivalent using a seasonal strength measure. The answer of one each is exactly the airline model specification, ARIMA(0,1,1)(0,1,1) with period 12, fitted in Exercise 5.7. Order matters in practice: take the seasonal difference first, then re-test, because removing seasonality sometimes removes the apparent trend as well and saves you an unnecessary ordinary difference. Over-differencing is a real cost, inflating variance and forcing spurious moving-average terms into the model.

Explanation: Two differences are required, not the one that most people assume from the plot, which is the practical lesson here. The counter cap is not decoration: without it, a series that never rejects sends the loop differencing until the vector is empty and R throws an obscure error deep inside the test. The 0.01 p-value is a floor, not an exact figure, because adf.test() interpolates from a printed table and warns when the statistic falls outside it. In production ndiffs() is preferable since it applies the same logic with better defaults, but writing the loop once makes clear what the helper is actually doing.

Explanation: KPSS reverses ADF: here the null is stationarity and a small p-value is evidence against it. At 0.054 the test sits right on the 5 percent line, so it barely fails to reject, while ADF failed to reject non-stationarity. Both point the same direction, towards a series that is not comfortably stationary, but neither is emphatic, which is normal for 100 observations. Use null = "Trend" when you suspect the series is stationary around a deterministic trend rather than a fixed level, since the two variants can easily disagree on a trending series.

Explanation: A p-value of 0.72 means there is no evidence of autocorrelation, which is the correct answer since the data really were drawn independently. Running a test against data whose truth you already know is the fastest way to calibrate your reading of it: compare 0.72 here with 1.255e-14 for Nile and the contrast is unmistakable. set.seed() is what makes the numbers reproducible, and without it every run gives a different statistic. This white-noise pattern is also the target for model residuals: a good model leaves behind something that looks exactly like this.

Explanation: Alpha controls how fast the level adapts. At 0.247 the model puts about a quarter of the weight on the newest observation and spreads the rest geometrically backwards, so the level responds to change but is not whipsawed by a single unusual flood year. Alpha near 1 would reduce the model to a naive last-value forecast, and near 0 would freeze it at the long-run mean. The single coefficient a of 805 is the final level estimate, and because there is no trend term, every future forecast from this model is that same flat 805.

Explanation: Both alpha and beta optimise to exactly 1, which is worth pausing on: the model has decided to put all weight on the most recent observation and the most recent slope, so it behaves like a random walk with drift. That happens when a series is dominated by a stochastic trend, and it is a signal to compare against an ARIMA fit rather than accept the smoothing model blindly. The b coefficient of -2 means the model expects a two-unit decline per period, so a 10-step forecast extrapolates 20 units below the current level of 220. Damping that extrapolation with forecast::holt(x, damped = TRUE) is usually safer over long horizons.

Explanation: Each parameter tells a different story. Alpha at 0.28 means the level updates moderately, beta at 0.03 means the trend is treated as almost fixed and changes only very slowly, and gamma at 0.87 means the seasonal factors are re-learned aggressively from each new year. That combination describes a series with a steady growth path and an evolving seasonal shape. seasonal = "multiplicative" is essential here, since the additive variant would force one fixed seasonal offset onto a swing that quadruples over the record. SSE is only comparable between models fitted to the same data on the same scale, so never compare it against a fit on log(AirPassengers).

Explanation: The three letters are error, trend and season, so ETS(M,Ad,M) means multiplicative errors, an additive damped trend, and multiplicative seasonality. The damping is why phi appears at 0.98, shrinking the extrapolated trend slightly at each step so long-horizon forecasts flatten instead of running away. Note gamma at 1e-04, effectively zero, meaning the seasonal shape is held fixed, which is the opposite of what Holt-Winters concluded in Exercise 5.3; the two routines optimise different likelihoods and often disagree. ets() selects by AICc, so its scores are comparable across ETS specifications but not against an ARIMA fitted to a transformed series.

Explanation: Judge each coefficient against roughly twice its standard error. The moving-average term at -0.8741 with a standard error of 0.0605 is unambiguously significant, while the autoregressive term at 0.2544 against 0.1194 is only about two standard errors out and is therefore borderline. That is a hint the simpler ARIMA(0,1,1) may fit nearly as well, which Exercise 5.8 confirms by comparing AIC across candidate orders. Because d = 1, no intercept is reported: differencing removes the mean, and R suppresses the constant unless you ask for drift.

Explanation: The search settles on one difference, not the two the ADF loop demanded, because auto.arima() picks d with a KPSS-based rule and then lets the ARMA terms absorb whatever non-stationarity is left. Disagreement between a mechanical differencing loop and the automatic selector is normal, and the AICc is what arbitrates. Two cautions. The default stepwise = TRUE explores only part of the order space, so auto.arima(x, stepwise = FALSE, approximation = FALSE) is worth running when accuracy matters more than speed. And models with different d are fitted to different data, so their AIC values are not comparable across differencing orders.

Explanation: This is the model Box and Jenkins built this dataset around, and it remains a strong baseline for seasonal series decades later. Each piece does one job: the log makes the widening seasonal swing additive, the ordinary difference removes the trend, the seasonal difference removes the annual pattern, and the two moving-average terms mop up what is left. Both coefficients are comfortably more than twice their standard errors, so neither is redundant. Note the sigma-squared of 0.001348 is on the log scale, which is why the AIC is negative and cannot be compared with a model fitted to raw counts.

Explanation: The top three are separated by under two AIC points, which is well inside the range where the differences are not meaningful, so the parsimonious ARIMA(0,1,1) is the sensible pick over the nominally best ARIMA(1,1,1). A gap of roughly 2 is negligible, 4 to 7 is moderate evidence, and above 10 is decisive. Note the row names survive the sort, which is why the index column runs 4, 5, 3, 2, 1 and doubles as a record of the original ordering. This comparison is only valid because every candidate uses d = 1 and is therefore fitted to identical data.

Explanation: This is regression with ARIMA errors: the step coefficient of -241.9, with a standard error of 33.0, estimates a drop of roughly 242 flow units after 1898 while the AR(1) term absorbs the remaining serial correlation. Modelling the break explicitly beats differencing it away, because you keep an interpretable estimate of the shift and you can quote its uncertainty. Compare the AIC of 1261.68 here with 1267.25 for ARIMA(1,1,1) in Exercise 5.5: the intervention model wins even after paying for an extra parameter. To forecast from this model you must supply future values of the regressor via newxreg, which for a permanent step is simply a run of 1s.

Explanation: The seasonal shape survives the round trip, with July and August forecast highest and November lowest, exactly as the historical pattern implies. The $mean element already carries the right calendar, so the exponentiated result prints as 1961 without any manual index work. One subtlety: exponentiating a forecast made on the log scale gives the median of the predicted distribution rather than the mean, so these numbers are slightly conservative as expectations. Exercise 6.9 shows the bias adjustment that converts them, and forecast(fit, h = 12, lambda = 0) applies the whole transform automatically.

Explanation: The point forecasts barely move, from 218.9 down to 217.0, while the interval width grows more than sevenfold from 12.4 to 92.9. That contrast is the real message of any forecast: the central path looks confident, but the honest uncertainty compounds with every step because each future error accumulates on top of the last. Reporting only the point forecast hides this entirely. Remember these intervals assume the model is correct and the residuals are normal, so in practice they are usually too narrow; forecast(f, bootstrap = TRUE) relaxes the normality assumption at least.

Explanation: Randomly sampling rows the way you would for cross-sectional data leaks the future into the training set and produces accuracy numbers that cannot be reproduced in production. A chronological split is the minimum defence. Two rules of thumb: the test set should be at least as long as the horizon you actually forecast, and for seasonal data it should span whole cycles, which is why 24 months rather than 20 is the right choice here. The lengths tying out to 120 plus 24 equals 144 confirms the split neither dropped nor duplicated a month.

Explanation: Read the diagnostics, not just the headline error. The mean error of 39.4 equals the mean absolute error exactly, which means every single forecast came in below the actual: the model is systematically under-predicting, not making random mistakes. The ACF1 of 0.46 confirms the errors are serially correlated rather than independent. Both symptoms trace back to the median-versus-mean issue from Exercise 6.1 plus genuine acceleration in 1959 to 1960 that the training window never saw. Theil's U below 1 still says the model beats a naive forecast, and RMSE punishes large misses harder than MAE, which is why it exceeds MAE here.

Explanation: The naive method repeats the last observed value forever and the seasonal naive repeats the value from the same month last year, which is why the seasonal version nearly halves every error measure on a strongly seasonal series. These numbers give the airline model from Exercise 6.4, at RMSE 43.2, a meaningful reference point: it beats seasonal naive by about 44 percent, so the extra complexity earns its keep. Any model that cannot beat snaive() on seasonal data should be discarded, and quoting a MAPE with no baseline beside it tells a stakeholder nothing about whether the forecast is good.

Explanation: The result matches the MAPE column in Exercise 6.4 to three decimals, which is the point: knowing the formula means you can compute the measure on any two vectors without a forecast object, and you can spot when a reported figure has been computed differently. MAPE has two well-known weaknesses. It is undefined when an actual is zero, and it penalises over-prediction more heavily than under-prediction because the denominator is the actual, which makes it a poor choice for intermittent demand. MASE, available in the same accuracy() output when a training series is supplied, avoids both problems by scaling against a naive benchmark.

Explanation: Here failing to reject is the good outcome: a p-value of 0.56 says the residuals are indistinguishable from white noise, so the model has left no systematic pattern behind. The fitdf = 2 is what drops the degrees of freedom from 24 to 22, accounting for the ma1 and sma1 coefficients; omitting it inflates the degrees of freedom and makes the test too lenient. Note the contrast with Exercise 6.4, where the holdout errors were clearly autocorrelated. Clean in-sample residuals plus dirty out-of-sample errors is the classic signature of a structural change after the training window, not of a badly specified model.

Explanation: tsCV() walks an expanding window forward through the series: it fits on everything up to time t, forecasts t plus 1, records the error, then moves on. That yields 96 out-of-sample errors from 100 observations rather than the single holdout of Exercise 6.4, so the resulting RMSE is a far more stable estimate. The early NA values appear because the first few windows are too short to fit an ARIMA. The cost is compute, since the model is refitted at every origin, which is why tsCV() is slow on long series and is usually run with a wrapper that reuses fitted coefficients instead.

Explanation: The adjustment is tiny here, about 0.07 percent, because the residual variance on the log scale is only 0.00135; on a noisier series it can move the forecast by several percent and matters a great deal when the forecasts are summed into a budget. The correction applies to the point forecast alone: quantiles are preserved under a monotone transform, so exponentiating the bounds is already exact and adjusting them would be wrong. In everyday work forecast(fit, h = 6, lambda = 0, biasadj = TRUE) does the whole log-and-back round trip for you, and lambda = "auto" picks a Box-Cox power instead of forcing a log.