ARMA, ARIMA and SARIMA Models

Stationarity, Identification, Estimation, Diagnostics, Forecasting and Validation

Published

December 3, 2023

Modified

September 15, 2026

Purpose

This section fits univariate models to four series: the Consumer Price Index, real GDP, and the quarterly averages of the 3-month and 6-month Treasury bill spreads over the Federal Funds Rate. Each series follows the same sequence. The order of differencing is set from the augmented Dickey-Fuller test and the autocorrelation function, the autoregressive and moving-average orders are identified from the ACF and PACF of the stationary series, a grid of candidate orders is estimated and ranked by AIC, AICc and BIC, and the ranking is compared with the selection made by auto.arima (Hyndman and Khandakar 2008). The retained model is then checked for coefficient significance and residual autocorrelation, used to forecast, compared with four benchmark methods, and, for the two spreads, tested by rolling-origin cross-validation.

The four series stand in different relations to the study’s question. The two spreads are the Federal Funds Rate relationship measured directly: each is a Treasury yield minus the policy rate on the same day, so a model of the spread is a model of how far the short end of the curve sits from where the Federal Reserve has set it. CPI and real GDP are the two variables the Federal Reserve’s mandate is written over, and the models here establish what each looks like on its own, before the ARIMAX and VAR sections add the policy rate to them. Nothing in this section uses the policy rate as an input; that is the next section’s subject.

Two points apply throughout. Information criteria are comparable only among models fitted to the same differenced series, so a model with d = 1 cannot be ranked against one with d = 2 by AIC; where the two appear together the comparison is made on diagnostics and out-of-sample error instead. The benchmark tables report training-set errors, which measure fit, not forecast accuracy; the cross-validation tables measure forecast accuracy.

Consumer Price Index

Data

The series is CPIAUCSL, monthly and seasonally adjusted, from January 1947 to August 2023, 920 observations. It is modelled in logarithms so that differences are approximately percentage changes and the variance does not grow with the level. Inflation is one of the Federal Reserve’s two mandate variables and the one the policy rate is most often moved to control; this section describes its own dynamics, and the ARIMAX section tests whether the policy rate adds to them.

Stationarity

Level                            ADF statistic  -0.82   p-value 0.960

Log level                        ADF statistic  -1.42   p-value 0.822

First difference of log          ADF statistic  -4.30   p-value < 0.01

Second difference of log         ADF statistic -15.04   p-value < 0.01

The level and the log level do not reject the unit root. The first difference of the log, which is monthly inflation, rejects it with a statistic of −4.30, but its ACF is significant at every one of the first 40 lags, with a lag-one value of 0.58 and a slow decline, and its PACF is significant through lag 7. Inflation is stationary by the test but so persistent that an ARMA model in first differences would need a high order to absorb the autocorrelation. The second difference has an ACF that cuts off after lag 2 (−0.31 and −0.16), apart from a cluster of small significant values at lags 10 to 12, and a PACF that tails off through lag 9. That pattern identifies a moving-average process of order 2. The model is therefore fitted with d = 2, which is also the order auto.arima selects.

Model Selection

The grid covers p from 0 to 3 and q from 0 to 2 with d = 2, twelve models, all without a constant since a constant with d = 2 would impose a quadratic trend.

Six Lowest-AICc Models
p	d	q	AIC	AICc	BIC
3	2	2	-8246.31	-8246.22	-8217.38
0	2	2	-8242.22	-8242.19	-8227.75
1	2	1	-8241.53	-8241.51	-8227.07
1	2	2	-8241.21	-8241.16	-8221.92
2	2	1	-8241.02	-8240.98	-8221.73
2	2	2	-8239.42	-8239.35	-8215.31

auto.arima selection:

Series: lcpi 
ARIMA(1,2,1)(2,0,0)[12] 

Coefficients:
         ar1      ma1     sar1     sar2
      0.2957  -0.8775  -0.2088  -0.1736
s.e.  0.0413   0.0217   0.0347   0.0362

sigma^2 = 6.977e-06:  log likelihood = 4149.36
AIC=-8288.72   AICc=-8288.66   BIC=-8264.61

AIC and AICc select ARIMA(3,2,2), at −8246.3; BIC selects ARIMA(0,2,2), at −8242.2 against −8217.4 for the larger model. The difference in AIC between the two is 4.1 points for three additional parameters. auto.arima, which also searches seasonal terms, selects SARIMA(1,2,1)(2,0,0)[12] with an AIC of −8288.7, lower than either. The seasonal autoregressive terms respond to the autocorrelation at lags 10 to 12 noted above, which remains in the series even though it is published seasonally adjusted.

Diagnostics

AIC -8246.31   AICc -8246.22   BIC -8217.38   residual s.d. 0.0027
Ljung-Box p-value at lag 8: 0.482   Ljung-Box p-value at lag 12: 0.000   Ljung-Box p-value at lag 20: 0.000

Coefficient	Estimate	Std. Error	z	p
ar1	1.1965	0.0650	18.42	0.0000
ar2	-0.3280	0.0541	-6.06	0.0000
ar3	0.0728	0.0367	1.98	0.0474
ma1	-1.7608	0.0567	-31.03	0.0000
ma2	0.7647	0.0545	14.04	0.0000

AIC -8242.22   AICc -8242.19   BIC -8227.75   residual s.d. 0.0027
Ljung-Box p-value at lag 8: 0.848   Ljung-Box p-value at lag 12: 0.000   Ljung-Box p-value at lag 20: 0.000

Coefficient	Estimate	Std. Error	z	p
ma1	-0.5605	0.0317	-17.70	0.0000
ma2	-0.2529	0.0330	-7.67	0.0000

AIC -8288.72   AICc -8288.66   BIC -8264.61   residual s.d. 0.0026
Ljung-Box p-value at lag 8: 0.678   Ljung-Box p-value at lag 12: 0.033   Ljung-Box p-value at lag 20: 0.014

Coefficient	Estimate	Std. Error	z
ar1	0.2957	0.0413	7.16
ma1	-0.8775	0.0217	-40.44
sar1	-0.2088	0.0347	-6.03
sar2	-0.1736	0.0362	-4.79

All coefficients in all three models are significant at the 5% level; the weakest is the third autoregressive term of ARIMA(3,2,2), with p = 0.047. The residual autocorrelation tells a consistent story. At lag 8 the Ljung-Box test does not reject white noise for any model. At lags 12 and 20 it rejects for both non-seasonal models, with p-values below 0.001, because of the autocorrelation at lags 10 to 12 that neither model addresses. The seasonal model reduces the problem, with p-values of 0.033 at lag 12 and 0.014 at lag 20, but does not remove it.

ARIMA(0,2,2) is retained for forecasting. It is the BIC choice, it has two parameters against four or five, and its short-lag residuals are clean. The remaining autocorrelation at lags 10 to 12 is small in magnitude and affects the forecast intervals more than the point forecast. The model is

\[(1 - B)^2 \log X_t = (1 + \theta_1 B + \theta_2 B^2)\, w_t, \qquad \hat\theta_1 = -0.56,\ \hat\theta_2 = -0.25,\]

which is equivalent to Holt’s linear-trend method: the forecast extends the level and slope estimated at the end of the sample.

Forecast

The forecast runs 60 months, from September 2023 to August 2028. The point forecast takes the index from 306.3 to 380.8, an average of 4.4% a year, above the 3.7% recorded over the twelve months to August 2023 because the slope is estimated from the 2021–2023 period. The 80% interval at the end of the horizon runs from 311 to 466 and the 95% interval from 280 to 518, in index points.

Benchmarks

Method	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
ARIMA(0,2,2)	0.0000	0.0027	0.0019	-0.0006	0.0450	0.0528	0.0009
Mean	0.0000	0.8627	0.7913	-4.1295	19.0820	22.4620	0.9974
Naive	0.0029	0.0045	0.0034	0.0678	0.0804	0.0966	0.5750
Seasonal naive	0.0344	0.0441	0.0352	0.7988	0.8224	1.0000	0.9849
Drift	0.0000	0.0034	0.0024	0.0002	0.0582	0.0684	0.5750

The errors are in log units on the training set. ARIMA(0,2,2) has the lowest RMSE, 0.0027, followed by the drift method at 0.0034, the naive method at 0.0045, the seasonal naive at 0.044 and the mean at 0.86. The ordering is the same on MAE, MAPE and MASE. The mean method has the smallest mean percentage error only because that measure is signed and its errors cancel.

Real GDP

Data

The series is GDPC1, quarterly, seasonally adjusted, from the first quarter of 1947 to the second quarter of 2023, 306 observations. It is modelled in logarithms. Output is the channel through which a change in the policy rate is expected to reach prices and employment, with a lag of several quarters; the macroeconomic vector autoregression in the next section estimates that response, and the model here provides the univariate baseline it improves on.

Stationarity

Level                            ADF statistic  -1.05   p-value 0.929

Log level                        ADF statistic  -1.12   p-value 0.918

First difference of log          ADF statistic  -7.36   p-value < 0.01

Second difference of log         ADF statistic -11.30   p-value < 0.01

The first difference of log GDP, which is quarterly real growth, rejects the unit root with a statistic of −7.36. Its ACF is significant at lags 1 and 2 only, with a lag-one value of 0.13, and its PACF at lag 1 only. The mean of the differenced series is 0.0076 per quarter, 3.1% a year, so the model needs a drift term. A second difference is not required: its ACF has a lag-one value of −0.49, the signature of a series that has been differenced once too often. The models are fitted with d = 1 and drift.

Model Selection

The grid covers p from 1 to 5 and q from 0 to 1 with d = 1 and d = 2, twenty models. Arima drops the drift term when d = 2, so those ten models have no constant.

Six Lowest-AICc Models
p	d	q	AIC	AICc	BIC
2	1	0	-1873.35	-1873.22	-1858.47
1	1	0	-1872.47	-1872.39	-1861.31
1	1	1	-1871.92	-1871.79	-1857.04
3	1	0	-1871.65	-1871.45	-1853.05
2	1	1	-1871.48	-1871.28	-1852.88
3	1	1	-1871.07	-1870.79	-1848.75

auto.arima selection:

Series: lgdp 
ARIMA(0,2,3)(1,0,1)[4] 

Coefficients:
          ma1      ma2      ma3     sar1    sma1
      -0.8748  -0.0010  -0.1093  -0.6082  0.5539
s.e.   0.0572   0.0756   0.0562   0.3597  0.3727

sigma^2 = 0.0001251:  log likelihood = 935.79
AIC=-1859.59   AICc=-1859.31   BIC=-1837.29

Every model with d = 1 ranks above every model with d = 2. AIC and AICc select ARIMA(2,1,0) with drift, at −1873.4; BIC selects ARIMA(1,1,0) with drift, at −1861.3, and that model is within one AIC point of the first. auto.arima selects SARIMA(0,2,3)(1,0,1)[4], with d = 2 and no drift; its AIC of −1859.6 is not comparable with the d = 1 models.

Diagnostics

AIC -1872.47   AICc -1872.39   BIC -1861.31   residual s.d. 0.0112
Ljung-Box p-value at lag 8: 0.607   Ljung-Box p-value at lag 12: 0.766   Ljung-Box p-value at lag 20: 0.858

Coefficient	Estimate	Std. Error	z	p
ar1	0.1336	0.0567	2.35	0.0186
drift	0.0076	0.0007	10.31	0.0000

AIC -1873.35   AICc -1873.22   BIC -1858.47   residual s.d. 0.0111
Ljung-Box p-value at lag 8: 0.838   Ljung-Box p-value at lag 12: 0.878   Ljung-Box p-value at lag 20: 0.934

Coefficient	Estimate	Std. Error	z	p
ar1	0.1208	0.0569	2.12	0.0339
ar2	0.0969	0.0569	1.70	0.0889
drift	0.0076	0.0008	9.36	0.0000

AIC -1859.59   AICc -1859.31   BIC -1837.29   residual s.d. 0.0112
Ljung-Box p-value at lag 8: 0.463   Ljung-Box p-value at lag 12: 0.748   Ljung-Box p-value at lag 20: 0.938

Coefficient	Estimate	Std. Error	z	p
ma1	-0.8748	0.0572	-15.28	0.0000
ma2	-0.0010	0.0756	-0.01	0.9891
ma3	-0.1093	0.0562	-1.95	0.0516
sar1	-0.6082	0.3597	-1.69	0.0909
sma1	0.5539	0.3727	1.49	0.1372

In ARIMA(1,1,0) both terms are significant: the autoregressive coefficient is 0.13 with p = 0.019 and the drift is 0.0076 with p below 0.001. The Ljung-Box p-values are 0.61, 0.77 and 0.86 at lags 8, 12 and 20. ARIMA(2,1,0) adds a second lag with p = 0.089 and does not improve the residuals. In the auto.arima model, four of the five coefficients have p-values above 0.05 and the second moving-average term is zero to three decimals, making the model over-parameterised. ARIMA(1,1,0) with drift is retained:

\[(1 - \phi B)(1 - B)\log X_t = c + w_t, \qquad \hat\phi = 0.13,\ \hat c = 0.0076.\]

Forecast

The horizon is 60 quarters, from the third quarter of 2023 to the second quarter of 2038. The point forecast grows at the drift rate, 3.1% a year, and reaches $35,057 billion in 2038. The 80% interval at that date runs from $30,855 to $39,832 billion and the 95% interval from $28,838 to $42,618. A forecast fifteen years out from a univariate model describes the average growth rate of the sample, not the path of the economy.

Benchmarks

Method	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
ARIMA(1,1,0)	0.0000	0.0111	0.0069	0.0014	0.0791	0.2021	-0.0149
Mean	0.0000	0.6857	0.5944	-0.6001	6.7265	17.2987	0.9903
Naive	0.0076	0.0136	0.0103	0.0864	0.1165	0.2995	0.1336
Seasonal naive	0.0306	0.0400	0.0344	0.3471	0.3890	1.0000	0.7982
Drift	0.0000	0.0112	0.0071	0.0012	0.0815	0.2080	0.1336

ARIMA(1,1,0) and the drift method are nearly identical, with RMSE of 0.0111 and 0.0112 and MAE of 0.0069 and 0.0071. They are the same model apart from the autoregressive correction. The naive method has RMSE 0.0136, the seasonal naive 0.040 and the mean 0.69.

3-Month Spread

Data

The series is T3MFF, the 3-month Treasury bill yield minus the effective Federal Funds Rate, daily from January 1982. The daily values are averaged within each calendar quarter, giving 168 quarters from the first quarter of 1982 to the fourth quarter of 2023. The final quarter contains only two trading days, 2 and 3 October 2023, and its value of 0.29 is the average of those two days. The quarterly series ranges from −1.49 points in the second quarter of 1982 to 0.55 in the third quarter of 2022, with a mean of −0.20. The mean is not stable across the sample: it is −0.62 over 1982–1989 and −0.10 from 1990 onward.

Stationarity

Level                            ADF statistic  -3.32   p-value 0.070

First difference                 ADF statistic  -7.29   p-value < 0.01

Level                            KPSS statistic   1.58   p-value < 0.01

Seasonal component range: -0.026 to 0.037 points. Series standard deviation: 0.333.

The evidence on the level is mixed at quarterly frequency. The ADF statistic is −3.32 with p = 0.070, a rejection at the 10% level but not at 5%; the KPSS test, whose null is stationarity, rejects with p below 0.01. On the daily data the ADF statistic is −8.28 and the unit root is rejected decisively (see the exploratory analysis). Averaging to quarters leaves 168 observations of a persistent series, and both tests have low power in that setting. The autocorrelation structure is unambiguous: the ACF declines geometrically, 0.82, 0.63, 0.47, 0.39, 0.35, and the PACF is significant at lag 1 only, at 0.82. That is the pattern of a stationary first-order autoregression, and a spread between two overnight-to-three-month rates is bounded by arbitrage, so it is modelled in levels.

Seasonality is negligible. The seasonal component from a classical decomposition ranges from −0.03 to +0.04 points against a standard deviation of 0.33, about a tenth of the variation. No seasonal differencing is applied; seasonal AR and MA terms are included in the grid so that the data can select them if they help.

Model Selection

The grid covers p and q from 0 to 3 and seasonal P and Q from 0 to 1, with d = D = 0 and a mean, 64 models.

Six Lowest-AICc Models
p	q	P	Q	AIC	AICc	BIC
0	3	1	1	-90.18	-89.48	-68.32
1	2	1	1	-89.93	-89.23	-68.06
1	0	0	0	-88.91	-88.77	-79.54
1	3	1	1	-89.52	-88.61	-64.52
1	3	0	0	-88.54	-88.02	-69.80
3	2	1	0	-88.77	-87.86	-63.78

Three Lowest-BIC Models
p	q	AIC	AICc	BIC
1	0	-88.91	-88.77	-79.54
2	0	-87.89	-87.65	-75.40
1	1	-87.75	-87.51	-75.26

auto.arima with d = 0:

Series: s3 
ARIMA(1,0,0)(2,0,1)[4] with non-zero mean 

Coefficients:
         ar1    sar1    sar2     sma1     mean
      0.8365  0.7918  0.1543  -0.8568  -0.2096
s.e.  0.0453  0.1059  0.0899   0.0728   0.1686

sigma^2 = 0.03253:  log likelihood = 50.64
AIC=-89.29   AICc=-88.77   BIC=-70.55

BIC selects AR(1), at −79.5, with the next model, AR(2), four points behind. AIC and AICc select ARMA(0,3)(1,1)[4] at −90.2 against −88.9 for AR(1), a gain of 1.3 points for four additional parameters. auto.arima constrained to d = 0 selects AR(1)(2,0,1)[4]; left unconstrained it selects SARIMA(3,1,0)(2,0,0)[4], differencing the series once, a choice the tests above do not support.

Diagnostics

AIC -88.91   AICc -88.77   BIC -79.54   residual s.d. 0.1828
Ljung-Box p-value at lag 8: 0.137   Ljung-Box p-value at lag 12: 0.152   Ljung-Box p-value at lag 20: 0.204

Coefficient	Estimate	Std. Error	z	p
ar1	0.8443	0.0421	20.06	0.0000
intercept	-0.2040	0.0873	-2.34	0.0194

AIC -89.29   AICc -88.77   BIC -70.55   residual s.d. 0.1804
Ljung-Box p-value at lag 8: 0.026   Ljung-Box p-value at lag 12: 0.063   Ljung-Box p-value at lag 20: 0.173

Coefficient	Estimate	Std. Error	z	p
ar1	0.8365	0.0453	18.49	0.0000
sar1	0.7918	0.1059	7.48	0.0000
sar2	0.1543	0.0899	1.72	0.0863
sma1	-0.8568	0.0728	-11.76	0.0000
intercept	-0.2096	0.1686	-1.24	0.2137

The AR(1) coefficient is 0.844 with a standard error of 0.042, three and a half standard errors below one, and the mean is −0.204 with p = 0.019. The residual standard deviation is 0.183 points. The residuals have one autocorrelation outside the confidence band, −0.19 at lag 3, which puts the Ljung-Box p-values at lags 3 and 4 at 0.024 and 0.035; from lag 5 onward they exceed 0.05, reaching 0.14, 0.15 and 0.20 at lags 8, 12 and 20. Models that add a third lag, AR(3) and ARMA(1,3), remove the lag-3 autocorrelation but their additional coefficients are not significant and both have higher AIC and BIC than AR(1). The seasonal alternative has a seasonal AR coefficient of 0.79 and a seasonal MA coefficient of −0.86, which nearly cancel, and its Ljung-Box p-value at lag 8 is 0.026, worse than the simpler model. AR(1) is retained:

\[X_t - \mu = \phi\,(X_{t-1} - \mu) + w_t, \qquad \hat\phi = 0.844,\ \hat\mu = -0.204.\]

The half-life of a shock is $\log 0.5 / \log 0.844 = 4.1$ quarters. The residuals are not normal (the largest, 0.81 points, is in the second quarter of 1982), so the prediction intervals, which assume normality, understate the tails.

Forecast

From the last value of 0.29 the point forecast falls to 0.21 in the first quarter of 2024, 0.05 by the fourth quarter of 2024 and −0.14 by the fourth quarter of 2026, converging on the sample mean. The 80% interval is −0.02 to 0.45 one quarter ahead and −0.57 to 0.29 twelve quarters ahead; the 95% interval at twelve quarters is −0.80 to 0.52. Given the shift in the mean after 1990, the convergence toward −0.20 rather than toward the post-1990 mean of −0.10 is a limitation of a constant-mean model.

Benchmarks

Method	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
AR(1)	0.0043	0.1817	0.1204	15.4035	112.5906	0.5525	0.0704
Mean	0.0000	0.3317	0.2505	3.2142	311.7334	1.1496	0.8236
Naive	0.0064	0.1888	0.1232	17.7233	112.8604	0.5653	-0.0293
Seasonal naive	0.0354	0.3283	0.2179	59.5800	334.7433	1.0000	0.6614
Drift	0.0000	0.1887	0.1234	20.7863	112.9409	0.5663	-0.0293

AR(1) has the lowest RMSE, 0.182, against 0.189 for the naive and drift methods, 0.328 for the seasonal naive and 0.332 for the mean. The percentage-error columns are not informative for this series because it crosses zero.

Cross-Validation

The rolling-origin comparison includes the retained AR(1), the seasonal alternative, and two models that difference the series, SARIMA(1,1,1)(2,1,1)[4] and the unconstrained auto.arima choice SARIMA(3,1,0)(2,0,0)[4]. The differenced models are included to test whether the decision not to difference costs forecast accuracy. The initial training set is the first 40% of the sample, 67 quarters to the third quarter of 1998, and each model is re-estimated at each of the 101 later quarters.

Initial training set: 67 quarters. One-step forecasts: 101. Four-step windows: 98.

One-Step-Ahead Errors
Model	MAE	MSE
AR(1)	0.1015	0.0197
AR(1)(2,0,1)[4]	0.0947	0.0191
SARIMA(1,1,1)(2,1,1)[4]	0.0901	0.0195
SARIMA(3,1,0)(2,0,0)[4]	0.0967	0.0226

Mean Absolute Error by Horizon
Model	h = 1	h = 2	h = 3	h = 4
AR(1)	0.1020	0.1565	0.1842	0.2118
AR(1)(2,0,1)[4]	0.0950	0.1418	0.1655	0.1897
SARIMA(1,1,1)(2,1,1)[4]	0.0916	0.1308	0.1437	0.1558
SARIMA(3,1,0)(2,0,0)[4]	0.0982	0.1408	0.1560	0.1680

One step ahead the four models are close: mean squared errors of 0.0197, 0.0191, 0.0195 and 0.0226, and mean absolute errors of 0.102, 0.095, 0.090 and 0.097. At longer horizons the differenced SARIMA(1,1,1)(2,1,1)[4] has the lowest error, 0.156 at four quarters against 0.212 for AR(1). The reason is the shift in the mean. The test period runs from late 1998, when the spread averaged −0.08, and includes 2009–2015, when it was within 0.1 of zero in 27 of 28 quarters; the AR(1), estimated on data that begin in 1982, pulls every forecast toward a mean that lies between −0.4 and −0.2 depending on the origin, while a differenced model carries the current level forward. The differenced model is not retained, because its first-order MA coefficient is estimated at −1.00, which is the model undoing the differencing, and its residuals fail the Ljung-Box test at lag 8. The comparison does show that the AR(1)’s constant mean is its weakest assumption for this series.

6-Month Spread

Data

The series is T6MFF, the 6-month bill yield minus the effective Federal Funds Rate, averaged by calendar quarter over the same 168 quarters, with the same two-day final quarter. It ranges from −1.03 points in the first quarter of 2008 to 1.05 in the third quarter of 2022, with a mean of −0.03: −0.26 over 1982–1989 and 0.03 from 1990.

Stationarity

Level                            ADF statistic  -3.81   p-value 0.020

First difference                 ADF statistic  -6.83   p-value < 0.01

Level                            KPSS statistic   0.61   p-value 0.021

Seasonal component range: -0.019 to 0.023 points. Series standard deviation: 0.332.

The ADF statistic is −3.81 with p = 0.020, rejecting the unit root at the 5% level; the KPSS test rejects stationarity with p = 0.021, so the two tests disagree at the same level. The ACF is 0.74, 0.55, 0.34, 0.21 and is not significant beyond lag 4; the PACF is significant at lag 1 only. The series decays faster than the 3-month spread. The seasonal component spans −0.02 to +0.02 points against a standard deviation of 0.33. The series is modelled in levels without seasonal differencing.

Model Selection

Six Lowest-AICc Models
p	q	P	Q	AIC	AICc	BIC
3	2	0	1	-27.82	-26.91	-2.82
3	2	1	0	-26.66	-25.76	-1.67
3	2	1	1	-26.06	-24.92	2.05
3	3	0	1	-25.88	-24.74	2.23
0	3	0	1	-24.38	-23.86	-5.64
3	2	0	0	-24.05	-23.35	-2.18

Three Lowest-BIC Models
p	P	Q	AIC	AICc	BIC
1	0	0	-21.59	-21.44	-12.21
1	0	1	-21.23	-20.99	-8.74
1	1	0	-20.77	-20.53	-8.28

auto.arima with d = 0:

Series: s6 
ARIMA(3,0,3)(1,0,1)[4] with zero mean 

Coefficients:
         ar1      ar2     ar3      ma1     ma2      ma3     sar1    sma1
      1.0636  -1.0257  0.5449  -0.3279  0.9266  -0.0701  -0.7204  0.7807
s.e.  0.1737   0.2735  0.2123   0.1896  0.1565   0.1373   0.3308  0.2749

sigma^2 = 0.04846:  log likelihood = 19.16
AIC=-20.32   AICc=-19.18   BIC=7.79

BIC selects AR(1), at −12.2. AIC and AICc select ARMA(3,2)(0,1)[4] at −27.8 against −21.6 for AR(1), with seven parameters against one. Among non-seasonal models AR(3) has AIC −22.5, one point below AR(1), and a worse BIC. auto.arima with d = 0 returns ARMA(3,3)(1,1)[4] with zero mean, eight parameters; unconstrained it returns SARIMA(0,1,0)(2,0,0)[4].

Diagnostics

AIC -21.59   AICc -21.44   BIC -12.21   residual s.d. 0.2237
Ljung-Box p-value at lag 8: 0.421   Ljung-Box p-value at lag 12: 0.615   Ljung-Box p-value at lag 20: 0.767

Coefficient	Estimate	Std. Error	z	p
ar1	0.7382	0.0514	14.37	0.0000
intercept	-0.0231	0.0645	-0.36	0.7206

AIC -22.53   AICc -22.16   BIC -6.91   residual s.d. 0.2217
Ljung-Box p-value at lag 8: 0.381   Ljung-Box p-value at lag 12: 0.582   Ljung-Box p-value at lag 20: 0.758

Coefficient	Estimate	Std. Error	z	p
ar1	0.7189	0.0775	9.28	0.0000
ar2	0.1645	0.1040	1.58	0.1136
ar3	-0.1846	0.0826	-2.23	0.0255
intercept	-0.0246	0.0556	-0.44	0.6576

The AR(1) coefficient is 0.738 with a standard error of 0.051. The mean is −0.023 with p = 0.72, not distinguishable from zero. The Ljung-Box p-values are 0.42, 0.62 and 0.77. In AR(3), the second lag has p = 0.11 and the third is −0.18 with p = 0.026; its residuals are no cleaner than the AR(1)’s. AR(1) is retained:

\[X_t - \mu = \phi\,(X_{t-1} - \mu) + w_t, \qquad \hat\phi = 0.738,\ \hat\mu = -0.023.\]

The half-life of a shock is 2.3 quarters, against 4.1 for the 3-month spread. The residual standard deviation is 0.224 points.

Forecast

From the last value of 0.25 the point forecast falls to 0.18 in the first quarter of 2024 and is within 0.02 of zero from the third quarter of 2025. The 80% interval twelve quarters ahead is −0.44 to 0.41 and the 95% interval is −0.67 to 0.63, essentially the unconditional distribution of the series.

Benchmarks

Method	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1
AR(1)	0.0005	0.2224	0.1515	75.0858	167.7207	0.5210	-0.0036
Mean	0.0000	0.3311	0.2327	161.2462	176.9307	0.8003	0.7392
Naive	0.0021	0.2388	0.1580	47.2249	225.3856	0.5433	-0.1434
Seasonal naive	0.0169	0.4100	0.2908	240.2875	374.3727	1.0000	0.6546
Drift	0.0000	0.2388	0.1581	42.2620	229.5025	0.5436	-0.1434

AR(1) has the lowest RMSE, 0.222, and MAE, 0.152; the naive and drift methods follow at 0.239 and 0.158, the mean at 0.331, and the seasonal naive at 0.410. The percentage-error columns are again uninformative.

Cross-Validation

The comparison includes AR(1), AR(3), and the two differenced models SARIMA(2,1,2)(2,1,1)[4] and SARIMA(0,1,0)(2,0,0)[4], the latter being the unconstrained auto.arima choice. The design is the same as for the 3-month spread.

Initial training set: 67 quarters. One-step forecasts: 101. Four-step windows: 98.

One-Step-Ahead Errors
Model	MAE	MSE
AR(1)	0.1159	0.0309
AR(3)	0.1197	0.0315
SARIMA(2,1,2)(2,1,1)[4]	0.1323	0.0350
SARIMA(0,1,0)(2,0,0)[4]	0.1180	0.0307

Mean Absolute Error by Horizon
Model	h = 1	h = 2	h = 3	h = 4
AR(1)	0.1174	0.1679	0.1925	0.2039
AR(3)	0.1199	0.1663	0.1933	0.2091
SARIMA(2,1,2)(2,1,1)[4]	0.1335	0.1778	0.2042	0.2112
SARIMA(0,1,0)(2,0,0)[4]	0.1199	0.1853	0.2308	0.2558

AR(1) has the lowest one-step mean absolute error, 0.116, against 0.120 for AR(3), 0.118 for SARIMA(0,1,0)(2,0,0)[4] and 0.132 for SARIMA(2,1,2)(2,1,1)[4]. At four quarters ahead the errors are 0.204 for AR(1), 0.209 for AR(3), 0.211 for SARIMA(2,1,2)(2,1,1)[4] and 0.256 for SARIMA(0,1,0)(2,0,0)[4]. AR(1) is the lowest or within 0.002 of the lowest at every horizon, and the differenced models fall behind as the horizon lengthens. For this series the mean is close to zero in both halves of the sample, so the constant-mean assumption that handicapped the 3-month model does not bind.

Summary

Retained Models
Series	Frequency	Differencing	Retained model	Selected by	Residual s.d.	Ljung-Box (lag 8, 12, 20)	Training RMSE vs. naive
CPI (log)	Monthly	d = 2	ARIMA(0,2,2)	BIC	0.0027	0.85, < 0.01, < 0.01	0.0027 vs. 0.0045
Real GDP (log)	Quarterly	d = 1, drift	ARIMA(1,1,0)	BIC	0.0112	0.61, 0.77, 0.86	0.0111 vs. 0.0136
3-month spread	Quarterly	None	AR(1)	BIC	0.183	0.14, 0.15, 0.20	0.182 vs. 0.189
6-month spread	Quarterly	None	AR(1)	BIC	0.224	0.42, 0.62, 0.77	0.222 vs. 0.239

Model Choice

The four series divide into two groups by their order of integration, and the model that fits each group follows from that. CPI and real GDP are trending series that fail the unit root test in levels; the spreads are bounded series that pass it. The two trending series are modelled in differences and the two spreads in levels, and in every case BIC chose the smallest model in the grid that the diagnostics could support.

For CPI the retained model is ARIMA(0,2,2) on the log index. The second difference is used because monthly inflation, the first difference, is stationary by the ADF test but autocorrelated at every lag out to 40; a model in first differences would have needed a high order to absorb that. The two moving-average coefficients, −0.56 and −0.25, make the model equivalent to Holt’s linear-trend method, so its forecast extends the level and slope estimated at the end of the sample. That slope reflects the inflation of 2021–2023 and produces a five-year path of 4.4% a year, above the 3.7% recorded over the twelve months to August 2023.

For real GDP the retained model is ARIMA(1,1,0) with drift on the log level. Quarterly real growth is stationary with a lag-one autocorrelation of 0.13 and no structure beyond lag 2, so a random walk with drift and one autoregressive correction is sufficient. The drift, 0.0076 per quarter, is 3.1% a year, and the forecast grows at that rate. This model replaced the ARIMA(1,1,1) that a grid restricted to q = 1 would have selected; the moving-average term in that model had p = 0.16 and its removal lowers every information criterion.

For both spreads the retained model is AR(1) with a mean. The ACF declines geometrically and the PACF cuts off at lag 1 in both series, and BIC prefers AR(1) over every larger model by more than three points. The autoregressive coefficients, 0.84 for the 3-month spread and 0.74 for the 6-month, give half-lives of 4.1 and 2.3 quarters: a shock to the 3-month spread takes about a year to decay by half, and a shock to the 6-month spread about seven months. The 6-month spread’s mean is not distinguishable from zero. Seasonal terms were offered to the grid and did not survive: where they enter the AIC-minimising models they appear as near-cancelling pairs, and the seasonal component from a classical decomposition is a tenth of the series’ standard deviation.

Residual Diagnostics

Three of the four retained models have residuals that pass the Ljung-Box test at every lag examined. The exception is CPI, whose residuals fail at lags 12 and 20 with p-values below 0.001 because of autocorrelation at lags 10 to 12. The series is published seasonally adjusted, but the adjustment leaves a small annual pattern in the second difference. auto.arima’s seasonal model, SARIMA(1,2,1)(2,0,0)[12], reduces the problem but does not remove it, and was not retained because the gain does not justify three additional parameters for a pattern this small. The 3-month spread’s AR(1) has one significant residual autocorrelation at lag 3, −0.19, which no model in the grid removes without adding insignificant coefficients. Residuals of both spread models are heavier-tailed than normal, with the largest values in 1982 and 2008, so the prediction intervals understate the tails.

Forecast Accuracy

Against the benchmark methods on the training set, the retained model has the lowest RMSE for all four series. The margins differ. For CPI the model’s RMSE is 60% of the naive method’s and 80% of the drift method’s. For GDP the model and the drift method are within 1% of each other, which is expected since the model is a drift model with a small correction. For the spreads the model’s RMSE is 4% and 7% below the naive method’s, a small margin, and the mean-percentage measures are not informative because both series cross zero.

The rolling-origin cross-validation on the spreads gives the only out-of-sample evidence. For the 6-month spread AR(1) has the lowest error at every horizon from one to four quarters, with a one-step MAE of 0.116 points. For the 3-month spread AR(1) is competitive one step ahead, with a mean squared error of 0.0197 against 0.0195 for the best alternative, but at four quarters ahead its MAE of 0.212 is above the 0.156 of a differenced SARIMA(1,1,1)(2,1,1)[4]. The cause is a shift in the spread’s mean, from −0.62 over 1982–1989 to −0.10 since 1990, which a constant-mean model cannot accommodate: every forecast is pulled toward a mean estimated from a sample that includes the 1980s. The differenced model avoids that by carrying the current level forward, but it does so through a moving-average coefficient of −1.00, which is the estimator undoing the differencing, and its residuals fail the Ljung-Box test at lag 8. The AR(1) is retained as the correctly identified model, with the caveat that its multi-quarter forecasts for the 3-month spread would be improved by a mean estimated on the post-1990 sample or by a model that allows the mean to shift.

The Federal Funds Rate

The policy rate does not appear as a variable in this section, but two of the four series are defined by it. The 3-month and 6-month spreads are yields net of the policy rate, so the AR(1) results say something about the relationship directly: a shock to the gap between the short end of the curve and the Federal Reserve’s rate has a half-life of 4.1 quarters at three months and 2.3 at six, and the gap has no trend, which is what arbitrage between an overnight rate and a bill of a few months’ maturity implies. The estimated mean of each spread is the average distance the market has kept from the policy rate over forty years, −0.20 points at three months and indistinguishable from zero at six. The shift in the 3-month mean, from −0.62 in the 1980s to −0.10 from 1990, is the clearest sign in this section that the relationship between market rates and the policy rate has not been constant.

CPI and real GDP enter as the mandate variables the policy rate is set to influence. Their univariate models are the baselines the following sections measure against when the policy rate is added.

Limitations

All four models are univariate. They describe how each series moves relative to its own past and say nothing about the Federal Funds Rate, which is the variable this project is about. The ARIMAX section that follows adds the policy rate as a regressor. Three further limits apply. The GDP forecast runs fifteen years, far beyond the horizon at which a univariate model is informative; it reports the sample’s average growth rate, not a projection of the economy. The CPI forecast inherits the 2021–2023 slope and should be read as a continuation of recent inflation rather than an expectation about it. And the final quarterly observation for both spreads is an average of two trading days, 2 and 3 October 2023, which adds noise to the last point the models fit and forecast from.

References

Hyndman, Rob J, and Yeasmin Khandakar. 2008. “Automatic Time Series Forecasting: The Forecast Package for R.” Journal of Statistical Software 27 (3): 1–22. https://doi.org/10.18637/jss.v027.i03.

--- title: "ARMA, ARIMA and SARIMA Models" subtitle: "Stationarity, Identification, Estimation, Diagnostics, Forecasting and Validation" date: "2023-12-3" date-modified: today date-format: long format: html: smooth-scroll: true toc: true code-fold: true code-tools: true lightbox: true embed-resources: true bibliography: citation.bib execute: echo: false warning: false message: false --- ```{r setup} suppressMessages({ library(readxl); library(dplyr); library(lubridate); library(ggplot2) library(forecast); library(tseries); library(astsa); library(knitr) }) knitr::opts_chunk$set(dev.args = list(bg = "#f4f4f4")) theme_set(theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e0e0e0", linewidth = 0.3), plot.background = element_rect(fill = "#f4f4f4", colour = NA), panel.background = element_rect(fill = NA, colour = NA), legend.background = element_rect(fill = NA, colour = NA), legend.key = element_rect(fill = NA, colour = NA), strip.background = element_rect(fill = NA, colour = NA), plot.title = element_text(size = 11, face = "plain"), axis.title = element_text(size = 9.5, colour = "#555555"), axis.text = element_text(size = 9, colour = "#555555"))) par(bg = "#f4f4f4") # Coefficient table, information criteria and Ljung-Box p-values for one fit. report <- function(fit, lags = c(8, 12, 20)) { se <- sqrt(diag(fit$var.coef)) z <- coef(fit) / se k <- length(coef(fit)) - any(names(coef(fit)) %in% c("intercept", "drift")) tab <- data.frame(Coefficient = names(coef(fit)), Estimate = round(coef(fit), 4), `Std. Error` = round(se, 4), z = round(z, 2), p = formatC(2 * pnorm(-abs(z)), format = "f", digits = 4), check.names = FALSE) lb <- sapply(lags, function(L) Box.test(residuals(fit), lag = L, type = "Ljung-Box", fitdf = k)$p.value) cat(sprintf("AIC %.2f AICc %.2f BIC %.2f residual s.d. %.4f\n", fit$aic, fit$aicc, fit$bic, sqrt(fit$sigma2))) cat(sprintf("Ljung-Box p-value at lag %d: %.3f", lags, lb), sep = " "); cat("\n") kable(tab, row.names = FALSE, align = "lrrrr") } # Residual diagnostic panel from astsa, without the optimiser trace. diag_plot <- function(x, p, d, q, P = 0, D = 0, Q = 0, S = -1) { invisible(capture.output(sarima(x, p, d, q, P, D, Q, S, details = TRUE))) } # Training-set accuracy of a fitted model against the four benchmark methods. bench <- function(fit, x, h, label, unit) { a <- rbind(accuracy(forecast(fit, h))[1, ], accuracy(meanf(x, h))[1, ], accuracy(naive(x, h))[1, ], accuracy(snaive(x, h))[1, ], accuracy(rwf(x, drift = TRUE, h))[1, ]) a <- data.frame(Method = c(label, "Mean", "Naive", "Seasonal naive", "Drift"), round(a, 4), check.names = FALSE) cat("\n\n"); print(kable(a, row.names = FALSE)); cat("\n\n") autoplot(x) + autolayer(meanf(x, h), series = "Mean", PI = FALSE) + autolayer(naive(x, h), series = "Naive", PI = FALSE) + autolayer(snaive(x, h), series = "Seasonal naive", PI = FALSE) + autolayer(rwf(x, h, drift = TRUE), series = "Drift", PI = FALSE) + autolayer(forecast(fit, h), series = label, PI = FALSE) + guides(colour = guide_legend(title = "Method")) + ggtitle("Benchmark Forecasts") + xlab("Year") + ylab(unit) } # Rolling-origin cross-validation. The first `frac` of the sample is the initial # training set; each later quarter is forecast 1 to `H` steps ahead from all data # before it, and the model is re-estimated at every origin. cv_compare <- function(x, specs, frac = 0.4, H = 4) { n <- length(x); split <- floor(n * frac); tsz <- n - split; st <- tsp(x)[1] + (split - 1) / frequency(x) one <- matrix(NA, tsz, length(specs)); multi <- lapply(specs, function(s) matrix(NA, tsz, H)) for (i in 1:tsz) { xtr <- window(x, end = st + (i - 1) / frequency(x)) xte1 <- x[split + i] xteH <- suppressWarnings(window(x, start = st + i / frequency(x), end = st + (i + H - 1) / frequency(x))) for (j in seq_along(specs)) { s <- specs[[j]] f <- tryCatch(Arima(xtr, order = s$order, seasonal = list(order = s$seasonal, period = frequency(x)), include.mean = s$mean, method = "ML"), error = function(e) NULL) if (is.null(f)) next fc <- forecast(f, h = H)$mean one[i, j] <- abs(fc[1] - xte1) if (length(xteH) == H) multi[[j]][i, ] <- abs(fc - xteH) } } nm <- sapply(specs, `[[`, "name") tab1 <- data.frame(Model = nm, MAE = round(colMeans(one, na.rm = TRUE), 4), MSE = round(colMeans(one^2, na.rm = TRUE), 4), check.names = FALSE) tabH <- data.frame(Model = nm, t(sapply(multi, function(m) round(colMeans(m, na.rm = TRUE), 4))), check.names = FALSE) names(tabH)[-1] <- paste0("h = ", 1:H) cat(sprintf("\n\nInitial training set: %d quarters. One-step forecasts: %d. Four-step windows: %d.\n\n", split, sum(!is.na(one[, 1])), sum(!is.na(multi[[1]][, 1])))) print(kable(tab1, row.names = FALSE, caption = "One-Step-Ahead Errors")); cat("\n\n") print(kable(tabH, row.names = FALSE, caption = "Mean Absolute Error by Horizon")); cat("\n\n") long <- do.call(rbind, lapply(seq_along(specs), function(j) data.frame(Model = nm[j], h = 1:H, MAE = colMeans(multi[[j]], na.rm = TRUE)))) ggplot(long, aes(h, MAE, colour = Model)) + geom_line() + geom_point() + scale_x_continuous(breaks = 1:H) + xlab("Forecast horizon (quarters)") + ylab("Mean absolute error") } adf_line <- function(x, label) { a <- suppressWarnings(adf.test(x)) cat(sprintf("%-32s ADF statistic %6.2f p-value %s\n", label, a$statistic, ifelse(a$p.value <= 0.01, "< 0.01", ifelse(a$p.value >= 0.99, "> 0.99", sprintf("%.3f", a$p.value))))) } ``` *** ### Purpose This section fits univariate models to four series: the Consumer Price Index, real GDP, and the quarterly averages of the 3-month and 6-month Treasury bill spreads over the Federal Funds Rate. Each series follows the same sequence. The order of differencing is set from the augmented Dickey-Fuller test and the autocorrelation function, the autoregressive and moving-average orders are identified from the ACF and PACF of the stationary series, a grid of candidate orders is estimated and ranked by AIC, AICc and BIC, and the ranking is compared with the selection made by `auto.arima` [@hyndman2008]. The retained model is then checked for coefficient significance and residual autocorrelation, used to forecast, compared with four benchmark methods, and, for the two spreads, tested by rolling-origin cross-validation. The four series stand in different relations to the study's question. The two spreads are the Federal Funds Rate relationship measured directly: each is a Treasury yield minus the policy rate on the same day, so a model of the spread is a model of how far the short end of the curve sits from where the Federal Reserve has set it. CPI and real GDP are the two variables the Federal Reserve's mandate is written over, and the models here establish what each looks like on its own, before the ARIMAX and VAR sections add the policy rate to them. Nothing in this section uses the policy rate as an input; that is the next section's subject. Two points apply throughout. Information criteria are comparable only among models fitted to the same differenced series, so a model with d = 1 cannot be ranked against one with d = 2 by AIC; where the two appear together the comparison is made on diagnostics and out-of-sample error instead. The benchmark tables report training-set errors, which measure fit, not forecast accuracy; the cross-validation tables measure forecast accuracy. ## Consumer Price Index ### Data The series is `CPIAUCSL`, monthly and seasonally adjusted, from January 1947 to August 2023, 920 observations. It is modelled in logarithms so that differences are approximately percentage changes and the variance does not grow with the level. Inflation is one of the Federal Reserve's two mandate variables and the one the policy rate is most often moved to control; this section describes its own dynamics, and the ARIMAX section tests whether the policy rate adds to them. ```{r} df <- na.omit(read_excel("data/data/cpi-u.xlsx")) cpi <- ts(as.numeric(df[[2]]), start = c(1947, 1), frequency = 12) lcpi <- log(cpi) autoplot(cpi) + ggtitle("CPI") + xlab("Year") + ylab("Index, 1982-84 = 100") ``` ### Stationarity ```{r} adf_line(cpi, "Level"); adf_line(lcpi, "Log level") adf_line(diff(lcpi), "First difference of log"); adf_line(diff(lcpi, differences = 2), "Second difference of log") ``` The level and the log level do not reject the unit root. The first difference of the log, which is monthly inflation, rejects it with a statistic of −4.30, but its ACF is significant at every one of the first 40 lags, with a lag-one value of 0.58 and a slow decline, and its PACF is significant through lag 7. Inflation is stationary by the test but so persistent that an ARMA model in first differences would need a high order to absorb the autocorrelation. The second difference has an ACF that cuts off after lag 2 (−0.31 and −0.16), apart from a cluster of small significant values at lags 10 to 12, and a PACF that tails off through lag 9. That pattern identifies a moving-average process of order 2. The model is therefore fitted with d = 2, which is also the order `auto.arima` selects. ```{r fig.width=9, fig.height=3} gridExtra::grid.arrange( ggAcf(lcpi, 40) + ggtitle("ACF, Log Level"), ggAcf(diff(lcpi), 40) + ggtitle("ACF, First Difference"), ggPacf(diff(lcpi), 40) + ggtitle("PACF, First Difference"), nrow = 1) gridExtra::grid.arrange( ggAcf(diff(lcpi, differences = 2), 40) + ggtitle("ACF, Second Difference"), ggPacf(diff(lcpi, differences = 2), 40) + ggtitle("PACF, Second Difference"), nrow = 1) ``` ### Model Selection The grid covers p from 0 to 3 and q from 0 to 2 with d = 2, twelve models, all without a constant since a constant with d = 2 would impose a quadratic trend. ```{r} grid <- expand.grid(p = 0:3, d = 2, q = 0:2) grid <- cbind(grid, t(apply(grid, 1, function(o) { m <- Arima(lcpi, order = o); c(AIC = m$aic, AICc = m$aicc, BIC = m$bic) }))) grid <- grid[order(grid$AICc), ]; grid[, 4:6] <- round(grid[, 4:6], 2) kable(head(grid, 6), row.names = FALSE, caption = "Six Lowest-AICc Models") cat("auto.arima selection:\n"); auto.arima(lcpi) ``` AIC and AICc select `ARIMA(3,2,2)`, at −8246.3; BIC selects `ARIMA(0,2,2)`, at −8242.2 against −8217.4 for the larger model. The difference in AIC between the two is 4.1 points for three additional parameters. `auto.arima`, which also searches seasonal terms, selects `SARIMA(1,2,1)(2,0,0)[12]` with an AIC of −8288.7, lower than either. The seasonal autoregressive terms respond to the autocorrelation at lags 10 to 12 noted above, which remains in the series even though it is published seasonally adjusted. ### Diagnostics ```{r fig.width=9, fig.height=6} diag_plot(lcpi, 3, 2, 2); report(Arima(lcpi, order = c(3, 2, 2))) ``` ```{r fig.width=9, fig.height=6} diag_plot(lcpi, 0, 2, 2); report(Arima(lcpi, order = c(0, 2, 2))) ``` ```{r fig.width=9, fig.height=6} diag_plot(lcpi, 1, 2, 1, 2, 0, 0, 12); report(Arima(lcpi, order = c(1, 2, 1), seasonal = c(2, 0, 0))) ``` All coefficients in all three models are significant at the 5% level; the weakest is the third autoregressive term of `ARIMA(3,2,2)`, with p = 0.047. The residual autocorrelation tells a consistent story. At lag 8 the Ljung-Box test does not reject white noise for any model. At lags 12 and 20 it rejects for both non-seasonal models, with p-values below 0.001, because of the autocorrelation at lags 10 to 12 that neither model addresses. The seasonal model reduces the problem, with p-values of 0.033 at lag 12 and 0.014 at lag 20, but does not remove it. `ARIMA(0,2,2)` is retained for forecasting. It is the BIC choice, it has two parameters against four or five, and its short-lag residuals are clean. The remaining autocorrelation at lags 10 to 12 is small in magnitude and affects the forecast intervals more than the point forecast. The model is $$(1 - B)^2 \log X_t = (1 + \theta_1 B + \theta_2 B^2)\, w_t, \qquad \hat\theta_1 = -0.56,\ \hat\theta_2 = -0.25,$$ which is equivalent to Holt's linear-trend method: the forecast extends the level and slope estimated at the end of the sample. ### Forecast ```{r} fit_cpi <- Arima(lcpi, order = c(0, 2, 2)) fc <- forecast(fit_cpi, h = 60) autoplot(fc) + ggtitle("Log CPI, 60-Month Forecast") + xlab("Year") + ylab("Log index") ``` The forecast runs 60 months, from September 2023 to August 2028. The point forecast takes the index from 306.3 to 380.8, an average of 4.4% a year, above the 3.7% recorded over the twelve months to August 2023 because the slope is estimated from the 2021–2023 period. The 80% interval at the end of the horizon runs from 311 to 466 and the 95% interval from 280 to 518, in index points. ### Benchmarks ```{r fig.width=9, fig.height=4} #| results: asis bench(fit_cpi, lcpi, 60, "ARIMA(0,2,2)", "Log index") ``` The errors are in log units on the training set. `ARIMA(0,2,2)` has the lowest RMSE, 0.0027, followed by the drift method at 0.0034, the naive method at 0.0045, the seasonal naive at 0.044 and the mean at 0.86. The ordering is the same on MAE, MAPE and MASE. The mean method has the smallest mean percentage error only because that measure is signed and its errors cancel. ## Real GDP ### Data The series is `GDPC1`, quarterly, seasonally adjusted, from the first quarter of 1947 to the second quarter of 2023, 306 observations. It is modelled in logarithms. Output is the channel through which a change in the policy rate is expected to reach prices and employment, with a lag of several quarters; the macroeconomic vector autoregression in the next section estimates that response, and the model here provides the univariate baseline it improves on. ```{r} df <- na.omit(read_excel("data/data/GDP.xlsx")) gdp <- ts(as.numeric(df[[2]]), start = c(1947, 1), frequency = 4) lgdp <- log(gdp) autoplot(gdp) + ggtitle("Real GDP") + xlab("Year") + ylab("$ billions") ``` ### Stationarity ```{r} adf_line(gdp, "Level"); adf_line(lgdp, "Log level") adf_line(diff(lgdp), "First difference of log"); adf_line(diff(lgdp, differences = 2), "Second difference of log") ``` The first difference of log GDP, which is quarterly real growth, rejects the unit root with a statistic of −7.36. Its ACF is significant at lags 1 and 2 only, with a lag-one value of 0.13, and its PACF at lag 1 only. The mean of the differenced series is 0.0076 per quarter, 3.1% a year, so the model needs a drift term. A second difference is not required: its ACF has a lag-one value of −0.49, the signature of a series that has been differenced once too often. The models are fitted with d = 1 and drift. ```{r fig.width=9, fig.height=3} gridExtra::grid.arrange( ggAcf(lgdp, 40) + ggtitle("ACF, Log Level"), ggAcf(diff(lgdp), 40) + ggtitle("ACF, First Difference"), ggPacf(diff(lgdp), 40) + ggtitle("PACF, First Difference"), nrow = 1) ``` ### Model Selection The grid covers p from 1 to 5 and q from 0 to 1 with d = 1 and d = 2, twenty models. `Arima` drops the drift term when d = 2, so those ten models have no constant. ```{r} grid <- expand.grid(p = 1:5, d = 1:2, q = 0:1) grid <- cbind(grid, t(apply(grid, 1, function(o) { m <- suppressWarnings(Arima(lgdp, order = o, include.drift = TRUE)); c(AIC = m$aic, AICc = m$aicc, BIC = m$bic) }))) grid <- grid[order(grid$AICc), ]; grid[, 4:6] <- round(grid[, 4:6], 2) kable(head(grid, 6), row.names = FALSE, caption = "Six Lowest-AICc Models") cat("auto.arima selection:\n"); auto.arima(lgdp) ``` Every model with d = 1 ranks above every model with d = 2. AIC and AICc select `ARIMA(2,1,0)` with drift, at −1873.4; BIC selects `ARIMA(1,1,0)` with drift, at −1861.3, and that model is within one AIC point of the first. `auto.arima` selects `SARIMA(0,2,3)(1,0,1)[4]`, with d = 2 and no drift; its AIC of −1859.6 is not comparable with the d = 1 models. ### Diagnostics ```{r fig.width=9, fig.height=6} diag_plot(lgdp, 1, 1, 0); report(Arima(lgdp, order = c(1, 1, 0), include.drift = TRUE)) ``` ```{r fig.width=9, fig.height=6} diag_plot(lgdp, 2, 1, 0); report(Arima(lgdp, order = c(2, 1, 0), include.drift = TRUE)) ``` ```{r fig.width=9, fig.height=6} diag_plot(lgdp, 0, 2, 3, 1, 0, 1, 4); report(Arima(lgdp, order = c(0, 2, 3), seasonal = c(1, 0, 1))) ``` In `ARIMA(1,1,0)` both terms are significant: the autoregressive coefficient is 0.13 with p = 0.019 and the drift is 0.0076 with p below 0.001. The Ljung-Box p-values are 0.61, 0.77 and 0.86 at lags 8, 12 and 20. `ARIMA(2,1,0)` adds a second lag with p = 0.089 and does not improve the residuals. In the `auto.arima` model, four of the five coefficients have p-values above 0.05 and the second moving-average term is zero to three decimals, making the model over-parameterised. `ARIMA(1,1,0)` with drift is retained: $$(1 - \phi B)(1 - B)\log X_t = c + w_t, \qquad \hat\phi = 0.13,\ \hat c = 0.0076.$$ ### Forecast ```{r} fit_gdp <- Arima(lgdp, order = c(1, 1, 0), include.drift = TRUE) fc <- forecast(fit_gdp, h = 60) autoplot(fc) + ggtitle("Log Real GDP, 60-Quarter Forecast") + xlab("Year") + ylab("Log $ billions") ``` The horizon is 60 quarters, from the third quarter of 2023 to the second quarter of 2038. The point forecast grows at the drift rate, 3.1% a year, and reaches \$35,057 billion in 2038. The 80% interval at that date runs from \$30,855 to \$39,832 billion and the 95% interval from \$28,838 to \$42,618. A forecast fifteen years out from a univariate model describes the average growth rate of the sample, not the path of the economy. ### Benchmarks ```{r fig.width=9, fig.height=4} #| results: asis bench(fit_gdp, lgdp, 60, "ARIMA(1,1,0)", "Log $ billions") ``` `ARIMA(1,1,0)` and the drift method are nearly identical, with RMSE of 0.0111 and 0.0112 and MAE of 0.0069 and 0.0071. They are the same model apart from the autoregressive correction. The naive method has RMSE 0.0136, the seasonal naive 0.040 and the mean 0.69. ## 3-Month Spread ### Data The series is `T3MFF`, the 3-month Treasury bill yield minus the effective Federal Funds Rate, daily from January 1982. The daily values are averaged within each calendar quarter, giving 168 quarters from the first quarter of 1982 to the fourth quarter of 2023. The final quarter contains only two trading days, 2 and 3 October 2023, and its value of 0.29 is the average of those two days. The quarterly series ranges from −1.49 points in the second quarter of 1982 to 0.55 in the third quarter of 2022, with a mean of −0.20. The mean is not stable across the sample: it is −0.62 over 1982–1989 and −0.10 from 1990 onward. ```{r} df <- na.omit(read_excel("data/data/3-mon-t-bill-rate.xlsx")) df$Date <- as.Date(df[[1]]); df$Value <- as.numeric(df[[2]]); df$Quarter <- floor_date(df$Date, "quarter") q <- df |> group_by(Quarter) |> summarize(Value = mean(Value), .groups = "drop") s3 <- ts(q$Value, start = c(1982, 1), frequency = 4) autoplot(s3) + ggtitle("3-Month Bill Minus Federal Funds Rate, Quarterly Average") + xlab("Year") + ylab("% points") ``` ### Stationarity ```{r} adf_line(s3, "Level"); adf_line(diff(s3), "First difference") k <- suppressWarnings(kpss.test(s3)); cat(sprintf("%-32s KPSS statistic %6.2f p-value %s\n", "Level", k$statistic, ifelse(k$p.value <= 0.01, "< 0.01", sprintf("%.3f", k$p.value)))) dc <- decompose(s3); cat(sprintf("Seasonal component range: %.3f to %.3f points. Series standard deviation: %.3f.\n", min(dc$figure), max(dc$figure), sd(s3))) ``` The evidence on the level is mixed at quarterly frequency. The ADF statistic is −3.32 with p = 0.070, a rejection at the 10% level but not at 5%; the KPSS test, whose null is stationarity, rejects with p below 0.01. On the daily data the ADF statistic is −8.28 and the unit root is rejected decisively (see the exploratory analysis). Averaging to quarters leaves 168 observations of a persistent series, and both tests have low power in that setting. The autocorrelation structure is unambiguous: the ACF declines geometrically, 0.82, 0.63, 0.47, 0.39, 0.35, and the PACF is significant at lag 1 only, at 0.82. That is the pattern of a stationary first-order autoregression, and a spread between two overnight-to-three-month rates is bounded by arbitrage, so it is modelled in levels. Seasonality is negligible. The seasonal component from a classical decomposition ranges from −0.03 to +0.04 points against a standard deviation of 0.33, about a tenth of the variation. No seasonal differencing is applied; seasonal AR and MA terms are included in the grid so that the data can select them if they help. ```{r fig.width=9, fig.height=3} gridExtra::grid.arrange( ggAcf(s3, 24) + ggtitle("ACF, Level"), ggPacf(s3, 24) + ggtitle("PACF, Level"), ggAcf(diff(s3), 24) + ggtitle("ACF, First Difference"), nrow = 1) ``` ### Model Selection The grid covers p and q from 0 to 3 and seasonal P and Q from 0 to 1, with d = D = 0 and a mean, 64 models. ```{r} sgrid <- function(x) { g <- expand.grid(p = 0:3, q = 0:3, P = 0:1, Q = 0:1) ic <- t(apply(g, 1, function(o) { m <- tryCatch(Arima(x, order = c(o["p"], 0, o["q"]), seasonal = c(o["P"], 0, o["Q"])), error = function(e) NULL) if (is.null(m)) c(AIC = NA, AICc = NA, BIC = NA) else c(AIC = m$aic, AICc = m$aicc, BIC = m$bic) })) g <- cbind(g, round(ic, 2)); g[order(g$AICc), ] } g3 <- sgrid(s3) kable(head(g3, 6), row.names = FALSE, caption = "Six Lowest-AICc Models") kable(head(g3[order(g3$BIC), ], 3), row.names = FALSE, caption = "Three Lowest-BIC Models") cat("auto.arima with d = 0:\n"); auto.arima(s3, d = 0) ``` BIC selects `AR(1)`, at −79.5, with the next model, `AR(2)`, four points behind. AIC and AICc select `ARMA(0,3)(1,1)[4]` at −90.2 against −88.9 for `AR(1)`, a gain of 1.3 points for four additional parameters. `auto.arima` constrained to d = 0 selects `AR(1)(2,0,1)[4]`; left unconstrained it selects `SARIMA(3,1,0)(2,0,0)[4]`, differencing the series once, a choice the tests above do not support. ### Diagnostics ```{r fig.width=9, fig.height=6} diag_plot(s3, 1, 0, 0); report(Arima(s3, order = c(1, 0, 0))) ``` ```{r fig.width=9, fig.height=6} diag_plot(s3, 1, 0, 0, 2, 0, 1, 4); report(Arima(s3, order = c(1, 0, 0), seasonal = c(2, 0, 1))) ``` The `AR(1)` coefficient is 0.844 with a standard error of 0.042, three and a half standard errors below one, and the mean is −0.204 with p = 0.019. The residual standard deviation is 0.183 points. The residuals have one autocorrelation outside the confidence band, −0.19 at lag 3, which puts the Ljung-Box p-values at lags 3 and 4 at 0.024 and 0.035; from lag 5 onward they exceed 0.05, reaching 0.14, 0.15 and 0.20 at lags 8, 12 and 20. Models that add a third lag, `AR(3)` and `ARMA(1,3)`, remove the lag-3 autocorrelation but their additional coefficients are not significant and both have higher AIC and BIC than `AR(1)`. The seasonal alternative has a seasonal AR coefficient of 0.79 and a seasonal MA coefficient of −0.86, which nearly cancel, and its Ljung-Box p-value at lag 8 is 0.026, worse than the simpler model. `AR(1)` is retained: $$X_t - \mu = \phi\,(X_{t-1} - \mu) + w_t, \qquad \hat\phi = 0.844,\ \hat\mu = -0.204.$$ The half-life of a shock is $\log 0.5 / \log 0.844 = 4.1$ quarters. The residuals are not normal (the largest, 0.81 points, is in the second quarter of 1982), so the prediction intervals, which assume normality, understate the tails. ### Forecast ```{r} fit_s3 <- Arima(s3, order = c(1, 0, 0)) fc <- forecast(fit_s3, h = 12) autoplot(fc) + ggtitle("3-Month Spread, 12-Quarter Forecast") + xlab("Year") + ylab("% points") ``` From the last value of 0.29 the point forecast falls to 0.21 in the first quarter of 2024, 0.05 by the fourth quarter of 2024 and −0.14 by the fourth quarter of 2026, converging on the sample mean. The 80% interval is −0.02 to 0.45 one quarter ahead and −0.57 to 0.29 twelve quarters ahead; the 95% interval at twelve quarters is −0.80 to 0.52. Given the shift in the mean after 1990, the convergence toward −0.20 rather than toward the post-1990 mean of −0.10 is a limitation of a constant-mean model. ### Benchmarks ```{r fig.width=9, fig.height=4} #| results: asis bench(fit_s3, s3, 12, "AR(1)", "% points") ``` `AR(1)` has the lowest RMSE, 0.182, against 0.189 for the naive and drift methods, 0.328 for the seasonal naive and 0.332 for the mean. The percentage-error columns are not informative for this series because it crosses zero. ### Cross-Validation The rolling-origin comparison includes the retained `AR(1)`, the seasonal alternative, and two models that difference the series, `SARIMA(1,1,1)(2,1,1)[4]` and the unconstrained `auto.arima` choice `SARIMA(3,1,0)(2,0,0)[4]`. The differenced models are included to test whether the decision not to difference costs forecast accuracy. The initial training set is the first 40% of the sample, 67 quarters to the third quarter of 1998, and each model is re-estimated at each of the 101 later quarters. ```{r fig.width=7, fig.height=3.5} #| results: asis cv_compare(s3, list( list(name = "AR(1)", order = c(1, 0, 0), seasonal = c(0, 0, 0), mean = TRUE), list(name = "AR(1)(2,0,1)[4]", order = c(1, 0, 0), seasonal = c(2, 0, 1), mean = TRUE), list(name = "SARIMA(1,1,1)(2,1,1)[4]", order = c(1, 1, 1), seasonal = c(2, 1, 1), mean = FALSE), list(name = "SARIMA(3,1,0)(2,0,0)[4]", order = c(3, 1, 0), seasonal = c(2, 0, 0), mean = FALSE))) ``` One step ahead the four models are close: mean squared errors of 0.0197, 0.0191, 0.0195 and 0.0226, and mean absolute errors of 0.102, 0.095, 0.090 and 0.097. At longer horizons the differenced `SARIMA(1,1,1)(2,1,1)[4]` has the lowest error, 0.156 at four quarters against 0.212 for `AR(1)`. The reason is the shift in the mean. The test period runs from late 1998, when the spread averaged −0.08, and includes 2009–2015, when it was within 0.1 of zero in 27 of 28 quarters; the `AR(1)`, estimated on data that begin in 1982, pulls every forecast toward a mean that lies between −0.4 and −0.2 depending on the origin, while a differenced model carries the current level forward. The differenced model is not retained, because its first-order MA coefficient is estimated at −1.00, which is the model undoing the differencing, and its residuals fail the Ljung-Box test at lag 8. The comparison does show that the `AR(1)`'s constant mean is its weakest assumption for this series. ## 6-Month Spread ### Data The series is `T6MFF`, the 6-month bill yield minus the effective Federal Funds Rate, averaged by calendar quarter over the same 168 quarters, with the same two-day final quarter. It ranges from −1.03 points in the first quarter of 2008 to 1.05 in the third quarter of 2022, with a mean of −0.03: −0.26 over 1982–1989 and 0.03 from 1990. ```{r} df <- na.omit(read_excel("data/data/6-mon-t-bill-rate.xlsx")) df$Date <- as.Date(df[[1]]); df$Value <- as.numeric(df[[2]]); df$Quarter <- floor_date(df$Date, "quarter") q <- df |> group_by(Quarter) |> summarize(Value = mean(Value), .groups = "drop") s6 <- ts(q$Value, start = c(1982, 1), frequency = 4) autoplot(s6) + ggtitle("6-Month Bill Minus Federal Funds Rate, Quarterly Average") + xlab("Year") + ylab("% points") ``` ### Stationarity ```{r} adf_line(s6, "Level"); adf_line(diff(s6), "First difference") k <- suppressWarnings(kpss.test(s6)); cat(sprintf("%-32s KPSS statistic %6.2f p-value %s\n", "Level", k$statistic, ifelse(k$p.value <= 0.01, "< 0.01", sprintf("%.3f", k$p.value)))) dc <- decompose(s6); cat(sprintf("Seasonal component range: %.3f to %.3f points. Series standard deviation: %.3f.\n", min(dc$figure), max(dc$figure), sd(s6))) ``` The ADF statistic is −3.81 with p = 0.020, rejecting the unit root at the 5% level; the KPSS test rejects stationarity with p = 0.021, so the two tests disagree at the same level. The ACF is 0.74, 0.55, 0.34, 0.21 and is not significant beyond lag 4; the PACF is significant at lag 1 only. The series decays faster than the 3-month spread. The seasonal component spans −0.02 to +0.02 points against a standard deviation of 0.33. The series is modelled in levels without seasonal differencing. ```{r fig.width=9, fig.height=3} gridExtra::grid.arrange( ggAcf(s6, 24) + ggtitle("ACF, Level"), ggPacf(s6, 24) + ggtitle("PACF, Level"), ggAcf(diff(s6), 24) + ggtitle("ACF, First Difference"), nrow = 1) ``` ### Model Selection ```{r} g6 <- sgrid(s6) kable(head(g6, 6), row.names = FALSE, caption = "Six Lowest-AICc Models") kable(head(g6[order(g6$BIC), ], 3), row.names = FALSE, caption = "Three Lowest-BIC Models") cat("auto.arima with d = 0:\n"); auto.arima(s6, d = 0) ``` BIC selects `AR(1)`, at −12.2. AIC and AICc select `ARMA(3,2)(0,1)[4]` at −27.8 against −21.6 for `AR(1)`, with seven parameters against one. Among non-seasonal models `AR(3)` has AIC −22.5, one point below `AR(1)`, and a worse BIC. `auto.arima` with d = 0 returns `ARMA(3,3)(1,1)[4]` with zero mean, eight parameters; unconstrained it returns `SARIMA(0,1,0)(2,0,0)[4]`. ### Diagnostics ```{r fig.width=9, fig.height=6} diag_plot(s6, 1, 0, 0); report(Arima(s6, order = c(1, 0, 0))) ``` ```{r fig.width=9, fig.height=6} diag_plot(s6, 3, 0, 0); report(Arima(s6, order = c(3, 0, 0))) ``` The `AR(1)` coefficient is 0.738 with a standard error of 0.051. The mean is −0.023 with p = 0.72, not distinguishable from zero. The Ljung-Box p-values are 0.42, 0.62 and 0.77. In `AR(3)`, the second lag has p = 0.11 and the third is −0.18 with p = 0.026; its residuals are no cleaner than the `AR(1)`'s. `AR(1)` is retained: $$X_t - \mu = \phi\,(X_{t-1} - \mu) + w_t, \qquad \hat\phi = 0.738,\ \hat\mu = -0.023.$$ The half-life of a shock is 2.3 quarters, against 4.1 for the 3-month spread. The residual standard deviation is 0.224 points. ### Forecast ```{r} fit_s6 <- Arima(s6, order = c(1, 0, 0)) fc <- forecast(fit_s6, h = 12) autoplot(fc) + ggtitle("6-Month Spread, 12-Quarter Forecast") + xlab("Year") + ylab("% points") ``` From the last value of 0.25 the point forecast falls to 0.18 in the first quarter of 2024 and is within 0.02 of zero from the third quarter of 2025. The 80% interval twelve quarters ahead is −0.44 to 0.41 and the 95% interval is −0.67 to 0.63, essentially the unconditional distribution of the series. ### Benchmarks ```{r fig.width=9, fig.height=4} #| results: asis bench(fit_s6, s6, 12, "AR(1)", "% points") ``` `AR(1)` has the lowest RMSE, 0.222, and MAE, 0.152; the naive and drift methods follow at 0.239 and 0.158, the mean at 0.331, and the seasonal naive at 0.410. The percentage-error columns are again uninformative. ### Cross-Validation The comparison includes `AR(1)`, `AR(3)`, and the two differenced models `SARIMA(2,1,2)(2,1,1)[4]` and `SARIMA(0,1,0)(2,0,0)[4]`, the latter being the unconstrained `auto.arima` choice. The design is the same as for the 3-month spread. ```{r fig.width=7, fig.height=3.5} #| results: asis cv_compare(s6, list( list(name = "AR(1)", order = c(1, 0, 0), seasonal = c(0, 0, 0), mean = TRUE), list(name = "AR(3)", order = c(3, 0, 0), seasonal = c(0, 0, 0), mean = TRUE), list(name = "SARIMA(2,1,2)(2,1,1)[4]", order = c(2, 1, 2), seasonal = c(2, 1, 1), mean = FALSE), list(name = "SARIMA(0,1,0)(2,0,0)[4]", order = c(0, 1, 0), seasonal = c(2, 0, 0), mean = FALSE))) ``` `AR(1)` has the lowest one-step mean absolute error, 0.116, against 0.120 for `AR(3)`, 0.118 for `SARIMA(0,1,0)(2,0,0)[4]` and 0.132 for `SARIMA(2,1,2)(2,1,1)[4]`. At four quarters ahead the errors are 0.204 for `AR(1)`, 0.209 for `AR(3)`, 0.211 for `SARIMA(2,1,2)(2,1,1)[4]` and 0.256 for `SARIMA(0,1,0)(2,0,0)[4]`. `AR(1)` is the lowest or within 0.002 of the lowest at every horizon, and the differenced models fall behind as the horizon lengthens. For this series the mean is close to zero in both halves of the sample, so the constant-mean assumption that handicapped the 3-month model does not bind. ## Summary ```{r} summary_tab <- data.frame( Series = c("CPI (log)", "Real GDP (log)", "3-month spread", "6-month spread"), Frequency = c("Monthly", "Quarterly", "Quarterly", "Quarterly"), `Differencing` = c("d = 2", "d = 1, drift", "None", "None"), `Retained model` = c("ARIMA(0,2,2)", "ARIMA(1,1,0)", "AR(1)", "AR(1)"), `Selected by` = c("BIC", "BIC", "BIC", "BIC"), `Residual s.d.` = c("0.0027", "0.0112", "0.183", "0.224"), `Ljung-Box (lag 8, 12, 20)` = c("0.85, < 0.01, < 0.01", "0.61, 0.77, 0.86", "0.14, 0.15, 0.20", "0.42, 0.62, 0.77"), `Training RMSE vs. naive` = c("0.0027 vs. 0.0045", "0.0111 vs. 0.0136", "0.182 vs. 0.189", "0.222 vs. 0.239"), check.names = FALSE) kable(summary_tab, row.names = FALSE, caption = "Retained Models") ``` ### Model Choice The four series divide into two groups by their order of integration, and the model that fits each group follows from that. CPI and real GDP are trending series that fail the unit root test in levels; the spreads are bounded series that pass it. The two trending series are modelled in differences and the two spreads in levels, and in every case BIC chose the smallest model in the grid that the diagnostics could support. For CPI the retained model is `ARIMA(0,2,2)` on the log index. The second difference is used because monthly inflation, the first difference, is stationary by the ADF test but autocorrelated at every lag out to 40; a model in first differences would have needed a high order to absorb that. The two moving-average coefficients, −0.56 and −0.25, make the model equivalent to Holt's linear-trend method, so its forecast extends the level and slope estimated at the end of the sample. That slope reflects the inflation of 2021–2023 and produces a five-year path of 4.4% a year, above the 3.7% recorded over the twelve months to August 2023. For real GDP the retained model is `ARIMA(1,1,0)` with drift on the log level. Quarterly real growth is stationary with a lag-one autocorrelation of 0.13 and no structure beyond lag 2, so a random walk with drift and one autoregressive correction is sufficient. The drift, 0.0076 per quarter, is 3.1% a year, and the forecast grows at that rate. This model replaced the `ARIMA(1,1,1)` that a grid restricted to q = 1 would have selected; the moving-average term in that model had p = 0.16 and its removal lowers every information criterion. For both spreads the retained model is `AR(1)` with a mean. The ACF declines geometrically and the PACF cuts off at lag 1 in both series, and BIC prefers `AR(1)` over every larger model by more than three points. The autoregressive coefficients, 0.84 for the 3-month spread and 0.74 for the 6-month, give half-lives of 4.1 and 2.3 quarters: a shock to the 3-month spread takes about a year to decay by half, and a shock to the 6-month spread about seven months. The 6-month spread's mean is not distinguishable from zero. Seasonal terms were offered to the grid and did not survive: where they enter the AIC-minimising models they appear as near-cancelling pairs, and the seasonal component from a classical decomposition is a tenth of the series' standard deviation. ### Residual Diagnostics Three of the four retained models have residuals that pass the Ljung-Box test at every lag examined. The exception is CPI, whose residuals fail at lags 12 and 20 with p-values below 0.001 because of autocorrelation at lags 10 to 12. The series is published seasonally adjusted, but the adjustment leaves a small annual pattern in the second difference. `auto.arima`'s seasonal model, `SARIMA(1,2,1)(2,0,0)[12]`, reduces the problem but does not remove it, and was not retained because the gain does not justify three additional parameters for a pattern this small. The 3-month spread's `AR(1)` has one significant residual autocorrelation at lag 3, −0.19, which no model in the grid removes without adding insignificant coefficients. Residuals of both spread models are heavier-tailed than normal, with the largest values in 1982 and 2008, so the prediction intervals understate the tails. ### Forecast Accuracy Against the benchmark methods on the training set, the retained model has the lowest RMSE for all four series. The margins differ. For CPI the model's RMSE is 60% of the naive method's and 80% of the drift method's. For GDP the model and the drift method are within 1% of each other, which is expected since the model is a drift model with a small correction. For the spreads the model's RMSE is 4% and 7% below the naive method's, a small margin, and the mean-percentage measures are not informative because both series cross zero. The rolling-origin cross-validation on the spreads gives the only out-of-sample evidence. For the 6-month spread `AR(1)` has the lowest error at every horizon from one to four quarters, with a one-step MAE of 0.116 points. For the 3-month spread `AR(1)` is competitive one step ahead, with a mean squared error of 0.0197 against 0.0195 for the best alternative, but at four quarters ahead its MAE of 0.212 is above the 0.156 of a differenced `SARIMA(1,1,1)(2,1,1)[4]`. The cause is a shift in the spread's mean, from −0.62 over 1982–1989 to −0.10 since 1990, which a constant-mean model cannot accommodate: every forecast is pulled toward a mean estimated from a sample that includes the 1980s. The differenced model avoids that by carrying the current level forward, but it does so through a moving-average coefficient of −1.00, which is the estimator undoing the differencing, and its residuals fail the Ljung-Box test at lag 8. The `AR(1)` is retained as the correctly identified model, with the caveat that its multi-quarter forecasts for the 3-month spread would be improved by a mean estimated on the post-1990 sample or by a model that allows the mean to shift. ### The Federal Funds Rate The policy rate does not appear as a variable in this section, but two of the four series are defined by it. The 3-month and 6-month spreads are yields net of the policy rate, so the `AR(1)` results say something about the relationship directly: a shock to the gap between the short end of the curve and the Federal Reserve's rate has a half-life of 4.1 quarters at three months and 2.3 at six, and the gap has no trend, which is what arbitrage between an overnight rate and a bill of a few months' maturity implies. The estimated mean of each spread is the average distance the market has kept from the policy rate over forty years, −0.20 points at three months and indistinguishable from zero at six. The shift in the 3-month mean, from −0.62 in the 1980s to −0.10 from 1990, is the clearest sign in this section that the relationship between market rates and the policy rate has not been constant. CPI and real GDP enter as the mandate variables the policy rate is set to influence. Their univariate models are the baselines the following sections measure against when the policy rate is added. ### Limitations All four models are univariate. They describe how each series moves relative to its own past and say nothing about the Federal Funds Rate, which is the variable this project is about. The ARIMAX section that follows adds the policy rate as a regressor. Three further limits apply. The GDP forecast runs fifteen years, far beyond the horizon at which a univariate model is informative; it reports the sample's average growth rate, not a projection of the economy. The CPI forecast inherits the 2021–2023 slope and should be read as a continuation of recent inflation rather than an expectation about it. And the final quarterly observation for both spreads is an average of two trading days, 2 and 3 October 2023, which adds noise to the last point the models fit and forecast from.