Time Series and ARIMA Method Card Template (Statistics, R and Python)

Time Series and ARIMA: [Your study question or course]

A method card for analysing and forecasting one time series with ARIMA and seasonal ARIMA (SARIMA): plot, decomposition, stationarity tests (ADF and KPSS), differencing, reading the ACF and PACF, model identification, fitting, residual checks, forecasting with intervals and honest accuracy measurement on a hold-out period. It is for students in a forecasting or econometrics course, thesis writers with monthly or weekly data, and analysts who must produce a forecast someone will plan around. It speeds up learning by giving you the Box and Jenkins steps in order, with a check at each step, and by running R and Python on the same series so you see where they differ by design. The worked examples use fictional data; replace them with your own. Source basis: Box and Jenkins (1970), Hyndman and Athanasopoulos (Forecasting: Principles and Practice), the NIST/SEMATECH e-Handbook (section 6.4) and the R forecast and statsmodels documentation.

Document control

When to use it

One-sentence rule: use ARIMA when you want to describe and forecast a single series from its own past, where consecutive observations are correlated and any trend or seasonality can be removed by differencing.

Ordinary regression assumes independent errors. In a monthly series, a busy month tends to be followed by another busy month, so that assumption fails and standard errors from ordinary regression are too small. ARIMA models the dependence directly: autoregressive (AR) terms use past values, moving average (MA) terms use past forecast errors, and the "I" (integrated) part is the number of differences needed to make the series stable. Seasonal ARIMA adds the same three pieces at the seasonal lag (12 for monthly data). Box and Jenkins (1970) turned this into a practical cycle of identification, estimation and diagnostic checking that is still the standard teaching route.

Alternatives: exponential smoothing (ETS), which often forecasts as well with fewer decisions; seasonal naive and drift benchmarks, which every model must beat; regression with ARIMA errors when you have predictors; STL plus a simple model for the seasonally adjusted series; and, for many related series, hierarchical or global models.

The question and the data

Question: how many visits will Harbour City Library receive each month in the next year, with a range the library manager can plan staffing around? A secondary question: is visiting growing, and how strong is the summer peak?

The series is the monthly count of entries recorded by the door counter at the main branch. It is a count, but with values in the thousands it behaves like a continuous variable. There are no gaps. Months 1 to 48 (2021 to 2024) are the training set; months 49 to 60 (2025) are the test set and are used only once, at the end.

The plot below shows the whole series, with the 2025 forecast from the model in Example A and its 95% prediction interval. Two features stand out before any calculation: a steady upward drift, and a strong yearly pattern with a July peak and a December low.

Assumptions and how to check each

How it works

Plain language first. A time series model says "this month is like recent months, plus a shock". The AR part says this month leans toward last month's value (or last year's same month). The MA part says that a surprise last month still echoes this month. Differencing removes slow movements: the first difference $y_t - y_{t-1}$ removes a wandering level, and the seasonal difference $y_t - y_{t-12}$ removes a stable yearly pattern. After differencing, what is left should look like a stable, patternless series with short memory, which a few AR and MA terms can describe.

Decomposition. STL (seasonal and trend decomposition using loess) splits the series into trend + seasonal + remainder. It is a description, not a forecast model, but it shows how strong each part is. A common summary is the strength of seasonality, $F_S = \max\!\left(0,\ 1 - \operatorname{Var}(R_t)/\operatorname{Var}(S_t + R_t)\right)$, and the same formula with the trend $T_t$ for trend strength. Values near 1 mean a strong component.

The model. With backshift operator $B$ (so $B y_t = y_{t-1}$), a seasonal ARIMA$(p,d,q)(P,D,Q)_m$ model is

where $\phi_p(B) = 1 - \phi_1 B - \dots - \phi_p B^p$ and $\theta_q(B) = 1 + \theta_1 B + \dots + \theta_q B^q$, and $\Phi$, $\Theta$ are the same polynomials in $B^m$. The model chosen in Example A, ARIMA$(1,0,0)(0,1,1)_{12}$ with drift, written out is

so the seasonally differenced series $w_t$ fluctuates around $12\delta$ (twelve months of drift $\delta$), with one month of AR memory and a seasonal MA term.

Stationarity tests. The augmented Dickey-Fuller (ADF) test has H0: the series has a unit root (needs differencing). KPSS reverses it, H0: the series is stationary (around a level or a trend). Use both: if ADF rejects and KPSS does not, the evidence agrees on "stationary". Both have low power with 48 points, so the plot and the seasonal strength matter as much as the p-values.

Reading ACF and PACF. The autocorrelation function (ACF) at lag $k$ is the correlation between $y_t$ and $y_{t-k}$. The partial autocorrelation (PACF) removes the effect of the lags in between. Approximate 95% bounds are $\pm 1.96/\sqrt{n}$. Rules of thumb on the differenced series: an AR($p$) shows a PACF that cuts off after lag $p$ and an ACF that decays; an MA($q$) shows an ACF that cuts off after lag $q$ and a PACF that decays; a spike at lag 12 in the ACF suggests a seasonal MA term, in the PACF a seasonal AR term. Real data rarely give a clean picture, so treat the rules as a way to propose two or three candidates, then compare them.

Accuracy. With forecast errors $e_t = y_t - \hat y_t$ on the test set of length $h$:

MAPE is undefined when any $y_t = 0$, explodes for small values, and penalises over- and under-forecasts unevenly; it is fine for this series (thousands of visits) but not for, say, daily counts of rare items. MASE below 1 means the model's test errors are smaller, on average, than the one-season-ahead errors of the seasonal naive method on the training data (the scale it is divided by).

Worked example A (full)

Goal: identify, fit and check a model on 2021 to 2024, then forecast 2025 and measure accuracy.

Step 1: decomposition. R code (install once with install.packages("forecast"); it brings urca):

visits <- c(4739, 5089, 5508, 5488, 5386, 5657, 5956, 5730, 5709, 5522, 5247, 4892,
            5259, 5285, 5427, 5702, 5647, 5865, 6296, 6023, 5605, 5723, 5404, 4970,
            5226, 5220, 5703, 5835, 5985, 6376, 6485, 6471, 6104, 5969, 5403, 5201,
            5476, 5544, 5853, 5934, 6001, 6269, 6525, 6461, 6356, 5972, 5416, 5223,
            5722, 5791, 5813, 5976, 6066, 6285, 6785, 6728, 6261, 6276, 5880, 5736)
y <- ts(visits, start = c(2021, 1), frequency = 12)
train <- window(y, end = c(2024, 12)); test <- window(y, start = c(2025, 1))
fit_stl <- stl(train, s.window = "periodic")
round(fit_stl$time.series[1:12, "seasonal"], 1)
rem <- fit_stl$time.series[, "remainder"]; tr <- fit_stl$time.series[, "trend"]
se <- fit_stl$time.series[, "seasonal"]
c(trend = 1 - var(rem) / var(tr + rem), season = 1 - var(rem) / var(se + rem))

Output (abridged and relabelled):

 [1] -449.1 -356.3  -34.9   70.4   73.7  352.3  617.7  465.2  229.3   72.5
[11] -366.2 -674.8
strength of trend: 0.817  seasonality: 0.937

The seasonal component (R, periodic, so the same every year) runs from about 675 visits below the trend in December to about 618 above it in July:

Seasonality strength 0.937 is very strong; trend strength 0.817 is strong. Python's STL(train, period=12) uses a seasonal smoother of length 7 by default, so its seasonal component can change from year to year and its first-year values differ (for example January -469.7 rather than -449.1). Neither is wrong: R's s.window = "periodic" forces a fixed pattern; set seasonal to a large odd number in Python to get close to that.

Step 2: stationarity and differencing.

library(urca); library(forecast)
summary(ur.df(train, type = "trend", lags = 2))       # ADF with trend
ur.kpss(train, type = "tau", lags = "short")          # KPSS around a trend
ndiffs(train); nsdiffs(train)
ds <- diff(train, lag = 12)                           # seasonal difference
ur.df(ds, type = "drift", lags = 1); ur.kpss(ds, type = "mu", lags = "short")
ndiffs(ds)                                            # after the seasonal difference

Output (abridged and relabelled):

ADF (trend, 2 lags) on train:  tau3 = -3.697   5% critical value -3.50
KPSS (trend) on train:         stat = 0.033    5% critical value 0.146
ndiffs: 1  nsdiffs: 1  ndiffs after seasonal diff: 0
ADF (drift, 1 lag) on seasonal difference:  tau2 = -3.325   5% critical value -2.93
KPSS (level) on seasonal difference:        stat = 0.225    5% critical value 0.463

Reading: on the raw series with a trend term, ADF rejects a unit root (-3.70 is below -3.50) and KPSS does not reject trend stationarity (0.033 is below 0.146), so a first difference is not needed if we allow a trend. nsdiffs says the seasonality is strong enough to take a seasonal difference ($D = 1$). (ndiffs says 1 only because by default it tests stationarity around a level, not a trend, which is exactly the trend in this series.) After seasonal differencing both tests agree the result is stationary, and ndiffs on it returns 0. Decision: $d = 0$, $D = 1$, with a drift term for the trend.

Step 3: ACF and PACF of the seasonally differenced series (36 values, bounds $\pm 1.96/\sqrt{36} = \pm 0.327$):

round(acf(ds, lag.max = 24, plot = FALSE)$acf[2:25], 2)
round(pacf(ds, lag.max = 24, plot = FALSE)$acf[1:24], 2)

Output (abridged and relabelled):

ACF  lags 1-12:  0.31  0.06  0.16  0.06  0.01  0.03 -0.05 -0.28 -0.19 -0.09 -0.12 -0.30
PACF lags 1-12:  0.31 -0.03  0.17 -0.05  0.01  0.01 -0.07 -0.27 -0.04 -0.03 -0.01 -0.28

Lag 1 is 0.31 in both, close to the bound, then both drop to near zero: this fits either AR(1) (PACF cuts off) or MA(1) (ACF cuts off), and with 36 values the two cannot be told apart. Lag 12 is negative in the ACF (-0.30), the classic sign of a seasonal MA(1) after seasonal differencing. Candidates: ARIMA$(1,0,0)(0,1,1)_{12}$ and ARIMA$(0,0,1)(0,1,1)_{12}$, both with drift.

Step 4: fit and check.

m1 <- Arima(train, order = c(1, 0, 0), seasonal = c(0, 1, 1), include.drift = TRUE)
m1
Box.test(residuals(m1), lag = 12, type = "Ljung-Box", fitdf = 2)
# Base R only, same model: drift as a regressor on time
m_base <- arima(train, order = c(1, 0, 0), seasonal = list(order = c(0, 1, 1), period = 12),
                xreg = 1:48)

Output (abridged and relabelled):

ARIMA(1,0,0)(0,1,1)[12] with drift
Coefficients:
         ar1     sma1    drift
      0.3611  -0.7360  14.7611
s.e.  0.1606   0.5612   1.9910
sigma^2 = 19068:  log likelihood = -231.12
AIC=470.23   AICc=471.53   BIC=476.57

Box-Ljung test: X-squared = 4.0217, df = 10, p-value = 0.9464

Reading: drift 14.76 visits per month (about 177 per year), SE 1.99, so growth is clearly positive. The AR(1) coefficient 0.36 means a busy month is followed by a slightly busy month. The seasonal MA coefficient -0.74 has a huge SE (0.56): with only four years, the seasonal term is poorly estimated. The Ljung-Box test on 12 lags, with df = 12 - 2 = 10 because two ARMA coefficients were estimated, gives p = .946: no evidence of leftover autocorrelation. Base R arima with xreg = 1:48 gives identical coefficients and log likelihood.

Step 5: forecast 2025 and measure accuracy.

f1 <- forecast(m1, h = 12)
accuracy(f1, test)
accuracy(snaive(train, h = 12), test)        # benchmark: same month last year

Output (abridged and rounded; the forecast table comes from forecast(m1, h = 12) printed next to test):

                 ME   RMSE    MAE   MPE MAPE MASE ACF1
Training set   4.24 114.50  77.02  0.03 1.34 0.40 0.00
Test set     -14.27 155.38 139.48 -0.22 2.32 0.73 0.32
Seasonal naive test set: RMSE 265.86, MAE 213.25, MAPE 3.52, MASE 1.11

         test  forecast  lo95  hi95
Jan 2025 5722      5576  5299  5854
Feb 2025 5791      5702  5408  5996
Mar 2025 5813      6053  5757  6349
...
Dec 2025 5736      5491  5195  5788

Python code (same data, same model):

import numpy as np, pandas as pd
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import adfuller, kpss, acf, pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox

visits = [4739, 5089, 5508, 5488, 5386, 5657, 5956, 5730, 5709, 5522, 5247, 4892,
          5259, 5285, 5427, 5702, 5647, 5865, 6296, 6023, 5605, 5723, 5404, 4970,
          5226, 5220, 5703, 5835, 5985, 6376, 6485, 6471, 6104, 5969, 5403, 5201,
          5476, 5544, 5853, 5934, 6001, 6269, 6525, 6461, 6356, 5972, 5416, 5223,
          5722, 5791, 5813, 5976, 6066, 6285, 6785, 6728, 6261, 6276, 5880, 5736]
y = pd.Series(visits, index=pd.period_range("2021-01", periods=60, freq="M"), dtype=float)
train, test = y[:48], y[48:]
print(adfuller(train, regression="ct", maxlag=2, autolag=None)[:2])
print(kpss(train, regression="ct", nlags="auto")[:2])
ds = train.diff(12).dropna()
print(np.round(acf(ds, nlags=12)[1:], 2), np.round(pacf(ds, nlags=12)[1:], 2))

t_train, t_test = np.arange(1, 49), np.arange(49, 61)       # drift as a regressor on time
m = SARIMAX(train, exog=t_train, order=(1, 0, 0), seasonal_order=(0, 1, 1, 12)).fit(disp=False)
print(m.summary().tables[1])
print(acorr_ljungbox(m.resid[12:], lags=[12], model_df=2))
fc = m.get_forecast(12, exog=t_test)
e = test - fc.predicted_mean
print(np.mean(np.abs(e)), np.sqrt(np.mean(e**2)), 100 * np.mean(np.abs(e / test)))
# Reproduce R exactly: likelihood of the differenced series
m_r = SARIMAX(train, exog=t_train, order=(1, 0, 0), seasonal_order=(0, 1, 1, 12),
              simple_differencing=True).fit(disp=False)
print(m_r.params, m_r.llf)

Output (abridged and relabelled; warnings omitted):

ADF (trend, 2 lags): stat = -3.697, p = 0.023
KPSS (trend): stat = 0.033, p = 0.100   (p-value is capped: actual p is larger)
ACF  1-12: [ 0.31  0.06  0.16  0.06  0.01  0.03 -0.05 -0.28 -0.19 -0.09 -0.12 -0.3 ]
PACF 1-12: [ 0.31 -0.04  0.19 -0.06  0.01  0.01 -0.09 -0.35 -0.06 -0.04 -0.03 -0.46]

                 coef    std err          z      P>|z|
x1            16.9704      3.289      5.159      0.000
ar.L1          0.3911      0.146      2.671      0.008
ma.S.L12      -0.3743      0.229     -1.632      0.103
sigma2        2.1e+04   6168.525      3.404      0.001
llf -232.08  AIC 472.15
Ljung-Box lag 12: stat = 3.389, p = 0.971

MAE = 150.65, RMSE = 175.99, MAPE = 2.51%

simple_differencing=True:  const 14.7613, ar.L1 0.3612, ma.S.L12 -0.7365, llf -231.12

Why R and Python differ here, by design. R's arima seasonally differences the series first and maximises the likelihood of the 36 differenced values. statsmodels SARIMAX by default keeps the differencing inside the state-space model and treats the first 12 values with a diffuse start, which is a different likelihood. With a poorly identified seasonal MA term, that small difference moves the estimate a lot (-0.74 in R, -0.37 in Python) and the drift with it (14.8 versus 17.0). Setting simple_differencing=True reproduces R to rounding (coefficients agree to about three decimals, log likelihood -231.12 in both). The PACF values also differ (most at lags 8 and 12) because the default estimators differ: R's pacf applies the Durbin-Levinson recursion to the usual (divide by n) ACF, while statsmodels pacf defaults to method="ywadjusted" (Yule-Walker with a divide by n - k ACF). Use method="ywm" in Python to get close to R. The forecasts are close in both cases (test MAE 139.5 in R, 150.7 in Python) and both beat the seasonal naive benchmark (MAE 213.3).

Calculation by hand for the key statistic. First seasonal difference: Jan 2022 minus Jan 2021 is $5259 - 4739 = 520$. Accuracy for the first three test months from the R forecasts: errors $5722 - 5576 = 146$, $5791 - 5702 = 89$, $5813 - 6053 = -240$. MAE over these three is $(146 + 89 + 240)/3 = 158.3$ and RMSE is $\sqrt{(146^2 + 89^2 + 240^2)/3} = \sqrt{28945.7} = 170.1$. RMSE is larger than MAE because squaring weights the March miss more heavily.

Effect size with CI. For a forecasting model, the "effect" is the drift and the size of the seasonal swing. Drift 14.76 visits per month, approximate 95% CI $14.76 \pm 1.96 \times 1.99 = [10.86, 18.66]$, so between about 130 and 224 extra visits per year. The seasonal swing from December to July is about 1,290 visits (STL). Forecast accuracy on 2025: MAE 139 visits, about 2.3% of a typical month (MAPE 2.32%), MASE 0.73.

Interpretation. Visits grew by roughly 180 a year, with a large and stable summer peak and a December low. After removing the yearly pattern, consecutive months are mildly correlated. The model forecast 2025 with a typical error of about 140 visits a month, 35% better than "same month last year" in MAE, and all 12 actual values fell inside the 95% prediction intervals (December 2025, at 5736, was closest to the edge, about 52 below the upper limit 5788; March 2025, at 5813, was about 56 above the lower limit 5757).

APA-style report sentence. A seasonal ARIMA(1,0,0)(0,1,1)[12] model with drift was fitted to monthly visits for 2021 to 2024 (n = 48; R forecast 9.0.2). Visits increased by 14.8 per month, 95% CI [10.9, 18.7]; residuals showed no remaining autocorrelation, Ljung-Box Q(10) = 4.02, p = .946. On the 2025 hold-out year the model achieved MAE = 139 visits, RMSE = 155, MAPE = 2.3% (MASE = 0.73), against MAE = 213 for the seasonal naive benchmark.

Worked example B (half solved)

auto.arima searches over orders and picks the model with the lowest AICc. Run it on the training set:

aa <- auto.arima(train)
aa
Box.test(residuals(aa), lag = 12, type = "Ljung-Box", fitdf = length(coef(aa)))
accuracy(forecast(aa, h = 12), test)

Output (abridged and relabelled):

ARIMA(0,0,1)(0,1,1)[12] with drift
Coefficients:
         ma1     sma1    drift
      0.4020  -0.7670  14.8232
s.e.  0.1584   0.6419   1.7866
sigma^2 = 18354:  log likelihood = -230.87
AIC=469.73   AICc=471.02   BIC=476.07

Box-Ljung test: X-squared = 3.9763, df = 9, p-value = 0.913
Test set: RMSE 153.48, MAE 135.10, MAPE 2.25, MASE 0.70

Python, for the same candidate: SARIMAX(train, exog=t_train, order=(0, 0, 1), seasonal_order=(0, 1, 1, 12)) gives AIC 472.09 (default likelihood), so it is a near tie with the AR(1) version in Python too.

1. Which non-seasonal term did auto.arima choose instead of AR(1)? [ ]
2. AICc of the manual model was 471.53 and of the automatic model [ ]. Is a difference of [ ] meaningful? [ ]
3. Both residual checks pass: manual p = .946, automatic p = [ ].
4. Test MAE: manual 139.48, automatic [ ]. Which would you report, and why does it barely matter? [ ]
5. Note one reason auto.arima can choose differently from a manual Box and Jenkins analysis: [ ] (fitdf in the Box.test call was [ ] because the drift is not an ARMA coefficient; strictly it should be 2 for both models, which changes the p-value only slightly.)

Exercise C (unsolved)

1. The seasonal naive forecast for each 2025 month is the same month in 2024. Using the data table, compute the 12 errors (actual 2025 minus 2024 value), then MAE and RMSE by hand or in a spreadsheet. Check with accuracy(snaive(train, h = 12), test) in R or a two-line numpy calculation.
2. In which months is the seasonal naive error largest, and why does a model with drift do better?
3. Compute MASE for the manual model, given that the in-sample seasonal naive MAE on the training set is 191.92 and the model's test MAE is 139.48.
4. Refit the manual model on all 60 months with Arima(y, order = c(1, 0, 0), seasonal = c(0, 1, 1), include.drift = TRUE) and look at the seasonal MA coefficient. What value is it close to, and what does that suggest about the seasonal pattern?

Reading software output

Reporting template

Fill-in sentence. "A seasonal ARIMA([p],[d],[q])([P],[D],[Q])[[m]] model [with drift / with a mean / without a constant] was fitted to [series, unit] for [training period] (n = [n]; [software and version]). [Main finding about trend or seasonality, with CI.] Residuals showed no remaining autocorrelation, Ljung-Box Q([df]) = [Q], p = [p]. On the hold-out period [test period] the model achieved MAE = [ ], RMSE = [ ], MAPE = [ ]% (MASE = [ ]), compared with [benchmark] (MAE = [ ])."

Results checklist.

• The series is plotted, with the train and test split marked.
• Transformations and differencing decisions are justified (plot, tests, seasonal strength).
• At least two candidate models were compared, plus the automatic choice.
• Residuals pass the Ljung-Box test with the correct df, and the residual ACF has no clear spikes.
• Accuracy is reported on the test set only, against a seasonal naive or naive benchmark.
• Final forecasts come from a model refitted on all data, with prediction intervals.
• Software, versions and any differences between R and Python are noted.

Common mistakes and misconceptions

• Random train and test split. Shuffling months into train and test lets the model see the future. Always split by time: train on the past, test on the most recent period.
• Reporting training-set accuracy. The training row of accuracy is always flattering. Report the test row.
• No benchmark. "MAPE 2.3%" means little until you know that seasonal naive gives 3.5%.
• Over-differencing. Taking both $d = 1$ and $D = 1$ when one is enough adds noise and pushes MA coefficients toward -1. A coefficient near -1 is a symptom.
• "Failing to reject KPSS proves stationarity." It does not; with 48 points both ADF and KPSS have low power. Use them alongside the plot.
• Comparing AIC across different differencing. AIC values are only comparable for models with the same $d$ and $D$, because the likelihood is computed on different data.
• Expecting R and Python to give the same coefficients. They maximise different likelihoods by default (see Example A). Check that forecasts, not every coefficient, agree.
• Ljung-Box with the wrong df. Subtract the number of estimated ARMA coefficients ($p + q + P + Q$), not the drift or mean.
• Treating a prediction interval as a confidence interval. A 95% prediction interval is for a single future month and includes the random shock; it is much wider than a CI for the mean forecast, and it ignores model uncertainty, so real coverage is often lower than 95%.
• Bootstrapping or regressing time series rows as if independent. See academic/stats-resampling-methods for the block bootstrap.

Self-test questions

1. The ACF of a differenced series has one spike at lag 1 and nothing after; the PACF decays slowly. Which model does this suggest?
2. ADF p = .40 and KPSS p < .01 on a monthly series. What do you conclude, and what do you do next?
3. Why must the test set come at the end of the series?
4. A model has MASE 1.15 on the test set. What does that mean for the library manager?
5. Your seasonal MA coefficient is -0.999. What is the likely problem, and what could you try?
6. Give one reason MAPE would be a poor accuracy measure for daily counts of interlibrary loan requests.

My notes and next steps

Sources and further reading

Structure aligned with the Box and Jenkins identification, estimation and diagnostic checking cycle; content based on the original method papers and the public-domain NIST/SEMATECH e-Handbook, written in our own words.

• Box, G. E. P., and Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day.
• Dickey, D. A., and Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366), 427–431.
• Kwiatkowski, D., Phillips, P. C. B., Schmidt, P., and Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54(1–3), 159–178.
• Ljung, G. M., and Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297–303.
• Cleveland, R. B., Cleveland, W. S., McRae, J. E., and Terpenning, I. (1990). STL: a seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6, 3–73.
• Hyndman, R. J., and Khandakar, Y. (2008). Automatic time series forecasting: the forecast package for R. Journal of Statistical Software, 27(3).
• Hyndman, R. J., and Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679–688.
• Hyndman, R. J., and Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed., online textbook, OTexts; further reading).
• NIST/SEMATECH e-Handbook of Statistical Methods, section 6.4 Introduction to Time Series Analysis, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• R forecast and urca documentation; statsmodels documentation (SARIMAX, STL, adfuller, kpss, acorr_ljungbox).
• Further reading: de Smith, M. J., Statistical Analysis Handbook (2024 edition, statsref.com); OpenIntro Statistics (CC BY-SA).

Answer key

Example B.

1. An MA(1) term (ma1 = 0.40) instead of AR(1).
2. 471.02; a difference of 0.51 AICc is negligible (differences below about 2 are not meaningful). The two models are practically equivalent: with $\phi = 0.36$, an AR(1) is close to an MA(1) with $\theta \approx 0.36$, since $1/(1 - 0.36B) \approx 1 + 0.36B + 0.13B^2 + \dots$.
3. p = .913.
4. 135.10. Either is fine; report the one chosen by your pre-stated rule (for example lowest AICc, so the automatic model), and note that the two forecasts differ by under 30 visits in every month (under 12 in most), far less than the forecast errors.
5. auto.arima uses a stepwise search, unit-root tests for $d$ and $D$, and AICc; the manual route reads the ACF and PACF, which cannot separate AR(1) from MA(1) with 36 values. The fitdf was 3 because length(coef(aa)) counted the drift.

Exercise C.

1. Errors (2025 minus 2024): 246, 247, -40, 42, 65, 16, 260, 267, -95, 304, 464, 513. Sum of absolute errors 2,559, MAE = 213.25. RMSE = 265.86 (R output).
2. November (464) and December (513), and in general most months are underpredicted: the seasonal naive method repeats last year, so it misses a whole year of growth (about 177 visits). A drift term adds that growth.
3. MASE = 139.48 / 191.92 = 0.73.
4. Refitted on 60 months: sma1 = -0.9987 (SE 1.06), drift 14.77, ar1 0.36. With $\Theta_1 \approx -1$ the MA factor $(1 - B^{12})$ on the right almost cancels the seasonal difference $(1 - B^{12})$ on the left, so the model is close to "fixed monthly means plus trend plus AR(1)": the seasonal pattern is essentially deterministic and the seasonal difference over-differenced the series. The estimate sits on the invertibility boundary (an MA root on the unit circle), where the usual standard error is unreliable (hence SE 1.06) and the model cannot be written as a convergent AR in past values. Try a regression on trend plus 11 monthly dummies with AR(1) errors (no seasonal difference), or ETS, and compare test accuracy. The R forecast for Jan 2026 from this refit is 5,898 visits, 95% PI [5,631, 6,165].

Self-test.

1. MA(1) on the differenced series: ARIMA(0,1,1) for the original series.
2. ADF cannot reject a unit root and KPSS rejects stationarity: both point to non-stationarity. Difference once (or seasonally difference if the problem is seasonal), then test again.
3. A forecast is used on future data; any test that lets the model see later values than the ones it forecasts overstates accuracy (look-ahead leakage).
4. The model's test MAE is 15% larger than the average error the "same month last year" method made on the training data (the MASE denominator). That is a warning sign that the model is no better than the simple method; compare it directly with the seasonal naive test MAE before telling the manager to switch.
5. The seasonal difference over-differenced a seasonal pattern that is essentially fixed (deterministic): the MA factor is cancelling the difference, and the model is at the non-invertible boundary, so its SE is unreliable. Drop the seasonal difference and use seasonal dummies (with ARMA errors) or seasonal AR terms; compare on the test set.
6. Many days have zero or very small counts, so $e_t / y_t$ is undefined or huge; MAPE would be dominated by a few small days. Use MAE or MASE.

Revision history