The fertility rate falls from about 125 births per 1,000 women in 1913 to the mid-60s by the early 1980s, with the post-war baby boom as a large hump in between. The sample autocorrelations decay extremely slowly:
| Lag \(k\) | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| \(\hat\rho_k\) | 0.945 | 0.873 | 0.807 | 0.734 | 0.649 | 0.551 |
Even six years apart, two observations still correlate above 0.5. This is a series with a great deal of memory: essentially all of the information in the recent past is still relevant today. A decay this gentle is also the pattern one would expect from a series with a unit root, or from one whose \(\rho\) is very close to one, which is what part (b) confirms.
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
fertil = pd.read_stata(
"https://empirical-economics.netlify.app/tutorials/datafiles/FERTIL3.DTA"
).sort_values("year")
fertil.plot(x="year", y="gfr")
plt.show()
plot_acf(fertil["gfr"], lags=6)
plt.show()| Quantity | Value |
|---|---|
| \(\hat\rho\) | 0.978 |
| Standard error | 0.026 |
| \(t\)-statistic | 37.60 |
| Half-life \(\ln(0.5)/\ln\hat\rho\) | 30.8 years |
| Implied long-run mean \(\hat\alpha/(1-\hat\rho)\) | 58.6 |
An AR(1) coefficient of 0.978 means that 98% of any deviation survives into the following year, and that half of a shock is still present after roughly 31 years. For an annual series with 72 observations, that is barely distinguishable from a shock that never dies at all.
So no, one should not be comfortable treating this series as stationary. The point estimate is less than one, but \((1-\hat\rho)/\text{s.e.}=(1-0.978)/0.026 \approx 0.86\), so the estimate sits less than one standard error away from one and the data cannot distinguish \(\rho=0.978\) from \(\rho=1\).
Treat that ratio as a back-of-the-envelope check rather than a test. A proper test of \(\rho=1\) is the Dickey–Fuller test, whose critical values are more negative than the normal ones precisely because the regressor is non-stationary under the null. Since we cannot reject even against the conventional threshold, we certainly cannot reject against the stricter one, so the verdict here is unchanged.
Note also that the implied long-run mean of 58.6 lies below every observation in the sample, which is a sign that reading a long-run mean off a near-unit-root AR(1) is not meaningful.
import numpy as np
import statsmodels.formula.api as smf
fertil["gfr_1"] = fertil["gfr"].shift(1)
ar1 = smf.ols("gfr ~ gfr_1", data=fertil).fit()
print(ar1.summary())
rho = ar1.params["gfr_1"]
print(np.log(0.5) / np.log(rho)) # half-life
print(ar1.params["Intercept"] / (1 - rho)) # implied long-run meanEvery model has to be estimated on the same observations, so all seven specifications drop the first six years and use the 66 remaining ones.
| Lags \(p\) | AIC | BIC |
|---|---|---|
| 0 | 576.95 | 581.33 |
| 1 | 387.08 | 393.65 |
| 2 | 382.55 | 391.31 |
| 3 | 382.46 | 393.40 |
| 4 | 378.15 | 391.29 |
| 5 | 377.17 | 392.50 |
| 6 | 374.86 | 392.38 |
AIC keeps falling all the way to \(p=6\) and would take more lags if we offered them. BIC reaches its minimum at \(p=4\), though \(p=2\) is only 0.02 higher, so BIC is effectively indifferent between two and four lags and clearly rejects six.
They disagree because they price a parameter differently. AIC charges \(2\) per parameter; BIC charges \(\ln n = \ln 66 \approx 4.19\), more than twice as much. Each extra lag has to improve the fit twice as much to be worth keeping under BIC, so BIC systematically selects the more parsimonious model. Neither answer is “correct”: AIC is aiming at forecast accuracy, BIC at identifying the true model if one is in the set considered.
One implementation detail to be aware of: the table above reports R’s values, and statsmodels will print numbers that are lower by a constant, because R counts the error variance as an estimated parameter and statsmodels does not. The constant is the same for every \(p\), so it cancels in any comparison and both languages select the same lag order. Never compare a criterion value across two different pieces of software.
max_p <- 6
lags <- sapply(1:max_p, function(k) dplyr::lag(fertil$gfr, k))
colnames(lags) <- paste0("L", 1:max_p)
# Drop the first max_p rows so that every model uses the same sample
d <- cbind(gfr = fertil$gfr, as.data.frame(lags))[-(1:max_p), ]
for (p in 0:max_p) {
f <- if (p == 0) gfr ~ 1 else
as.formula(paste("gfr ~", paste0("L", 1:p, collapse = " + ")))
m <- lm(f, data = d)
cat(p, round(AIC(m), 2), round(BIC(m), 2), "\n")
}from statsmodels.tsa.ar_model import ar_select_order
# ar_select_order handles the common-sample requirement internally
for ic in ["aic", "bic"]:
sel = ar_select_order(fertil["gfr"], maxlag=6, ic=ic, old_names=False)
print(ic, sel.ar_lags)
# The explicit version, for comparison with the R loop
max_p = 6
d = pd.DataFrame({"gfr": fertil["gfr"]})
for k in range(1, max_p + 1):
d[f"L{k}"] = fertil["gfr"].shift(k)
d = d.iloc[max_p:]
for p in range(max_p + 1):
rhs = " + ".join(f"L{k}" for k in range(1, p + 1)) or "1"
m = smf.ols(f"gfr ~ {rhs}", data=d).fit()
print(p, round(m.aic, 2), round(m.bic, 2))The AR(2) model estimated on the full sample is
\[ \widehat{gfr}_t = 2.636 + 1.271\,gfr_{t-1} - 0.305\,gfr_{t-2}. \]
The last two observations are \(gfr_{1984}=65.4\) and \(gfr_{1983}=65.8\). For 1985 both inputs are known:
\[ \hat g_{1985} = 2.636 + 1.271(65.4) - 0.305(65.8) = 65.69 . \]
For 1986 we do not know 1985, so we substitute the forecast we just computed, and 1984 moves into the second lag:
\[ \hat g_{1986} = 2.636 + 1.271(65.69) - 0.305(65.4) = 66.18 . \]
Note the direction: the forecast turns upward, even though the series has been falling. The sum of the AR coefficients is \(1.271-0.305=0.966\), so the model has an unconditional mean of \(2.636/(1-0.966)\approx 77.5\), well above the current level, and the forecast begins crawling towards it.
fertil <- fertil |> mutate(gfr_2 = lag(gfr, 2))
ar2 <- lm(gfr ~ gfr_1 + gfr_2, data = fertil)
b <- coef(ar2)
y_T <- fertil$gfr[nrow(fertil)] # 1984
y_T_1 <- fertil$gfr[nrow(fertil) - 1] # 1983
f1 <- b[1] + b[2] * y_T + b[3] * y_T_1 # 1985
f2 <- b[1] + b[2] * f1 + b[3] * y_T # 1986, using the 1985 forecast
c(f1, f2)
b[1] / (1 - b[2] - b[3]) # unconditional meanfrom statsmodels.tsa.ar_model import AutoReg
ar2 = AutoReg(fertil["gfr"], lags=2, old_names=False).fit()
b = ar2.params
y_T, y_T_1 = fertil["gfr"].iloc[-1], fertil["gfr"].iloc[-2]
f1 = b["const"] + b["gfr.L1"] * y_T + b["gfr.L2"] * y_T_1 # 1985
f2 = b["const"] + b["gfr.L1"] * f1 + b["gfr.L2"] * y_T # 1986
print(f1, f2)
print(b["const"] / (1 - b["gfr.L1"] - b["gfr.L2"])) # unconditional mean| Year | 1985 | 1986 | 1987 | 1988 | 1989 | 1990 | 1991 | 1992 | 1993 | 1994 |
|---|---|---|---|---|---|---|---|---|---|---|
| Forecast | 65.69 | 66.18 | 66.71 | 67.24 | 67.75 | 68.24 | 68.70 | 69.14 | 69.56 | 69.96 |
The forecast climbs steadily and is heading for the unconditional mean of about 77.5. What determines that destination is entirely the estimated coefficients: \(\hat\alpha/(1-\sum\hat\phi_i)\). Because \(\sum\hat\phi_i=0.966\) is so close to one, the approach is very slow: after ten years the forecast has covered under 40% of the distance from 65.4 to 77.5.
This is the general property of a stationary model: it can say something specific about the near future, but at long horizons it can only tell you the historical average. Here that average is not a credible statement about American fertility in 1994, which is a reason to distrust the model at long horizons rather than a reason to believe the forecast.
For a random walk \(Y_t=Y_{t-1}+u_t\), tomorrow’s value is \(Y_{t+1}=Y_t+u_{t+1}\). Taking the conditional expectation given everything known today,
\[ E[Y_{t+1}\mid I_t] = E[Y_t + u_{t+1}\mid I_t] = Y_t + E[u_{t+1}\mid I_t] = Y_t . \]
The single best piece of information is therefore today’s price, and nothing else in the history of the series adds anything. The whole past is summarised by the most recent value, which is why the random walk is the natural benchmark against which any financial forecasting model has to prove itself.
Iterating the process forward from \(t\),
\[ Y_{t+k} = Y_t + u_{t+1} + u_{t+2} + \dots + u_{t+k} . \]
The shock \(u_t\) is not in that sum, because it is already inside \(Y_t\). And \(Y_t\) enters \(Y_{t+k}\) with a coefficient of exactly one for every horizon \(k\), however large. So a piece of unexpected news shifts the entire future path of the price up or down by its full amount, permanently. There is no mechanism in the process that pulls the series back.
Contrast this with a stationary AR(1), where the same shock enters \(Y_{t+k}\) with coefficient \(\rho^{\,k}\), which goes to zero. The random walk is the case where that decay factor equals one.
Statistically the two are almost indistinguishable in this sample, which is exactly the point about the low power of unit-root tests. Economically they are different claims about the world.
If \(\rho=0.978\) the series has an equilibrium level. A shock is a temporary deviation, the half-life is 31 years, and a forecaster can in principle name a long-run value the series will return to. If \(\rho=1\) there is no equilibrium at all: the long-run mean and variance do not exist, every shock resets the level permanently, and the best long-horizon forecast is simply the last observation.
The practical consequence is that a policy question such as “will fertility recover on its own?” has opposite answers under the two models, while the data cannot decide between them. That is what you should report: the answer depends on an assumption the sample cannot verify, so state the assumption and check whether your conclusion survives both.
\[ \begin{aligned} E(Y_t) &= E(\mu + u_t + \theta u_{t-1}) \\ &= \mu + E(u_t) + \theta E(u_{t-1}) \\ &= \mu + 0 + 0 = \mu . \end{aligned} \]
A constant does not affect a variance, and \(u_t\) and \(u_{t-1}\) are uncorrelated because \(u\) is white noise, so the variance of the sum is the sum of the variances:
\[ \begin{aligned} \gamma_0 = \operatorname{Var}(Y_t) &= \operatorname{Var}(u_t + \theta u_{t-1}) \\ &= \operatorname{Var}(u_t) + \theta^2\operatorname{Var}(u_{t-1}) \\ &= \sigma^2_u + \theta^2\sigma^2_u = (1+\theta^2)\sigma^2_u . \end{aligned} \]
\[ \begin{aligned} \gamma_1 &= E\big[(Y_t-\mu)(Y_{t-1}-\mu)\big] \\ &= E\big[(u_t+\theta u_{t-1})(u_{t-1}+\theta u_{t-2})\big] \\ &= E[u_tu_{t-1}] + \theta E[u_tu_{t-2}] + \theta E[u_{t-1}^2] + \theta^2 E[u_{t-1}u_{t-2}] . \end{aligned} \]
White noise gives \(E[u_tu_s]=0\) for \(t\neq s\) and \(E[u_t^2]=\sigma^2_u\), so only the third term survives:
\[ \gamma_1 = \theta\sigma^2_u . \]
\[ \begin{aligned} \gamma_2 &= E\big[(u_t+\theta u_{t-1})(u_{t-2}+\theta u_{t-3})\big] \\ &= E[u_tu_{t-2}] + \theta E[u_tu_{t-3}] + \theta E[u_{t-1}u_{t-2}] + \theta^2 E[u_{t-1}u_{t-3}] . \end{aligned} \]
Every product now involves two different dates, so every expectation is zero:
\[ \gamma_2 = 0 . \]
The same argument applies to every \(k\ge 2\), so \(\gamma_k=0\) for all \(k\ge 2\).
An MA(1) process has a memory of exactly one period. A shock at \(t-2\) or earlier has no correlation whatsoever with \(Y_t\), and dividing through by \(\gamma_0\) gives
\[ \rho_1 = \frac{\theta}{1+\theta^2}, \qquad \rho_k = 0 \ \ (k\ge 2). \]
In a correlogram this looks completely different from an AR(1). The AR(1) shows \(\rho_k=\rho^{\,k}\), a smooth geometric decay that in principle never quite reaches zero. The MA(1) shows one spike at lag one and then a hard cut-off to zero. That contrast is the standard way of telling the two structures apart by eye.
| Quantity | Value |
|---|---|
Coefficient on pop |
0.872 |
| Standard error | 0.085 |
| \(t\)-statistic | 10.23 |
| \(R^2\) | 0.723 |
A naive reading is emphatic: population is a powerful driver of housing investment, explaining 72% of its variation with a \(t\)-statistic above ten. The coefficient says that a thousand extra people are associated with about 0.87 million 1982 dollars of extra housing investment, roughly $870 per additional person, which is well above average per-capita housing investment of about $520 in this sample. That discrepancy is a first hint that the regression is picking up something other than a marginal effect: both series rise steadily over the whole post-war period, and a regression of one trend on another will always look like this.
| Lags tested | LM statistic | \(p\)-value |
|---|---|---|
| 1 | 14.02 | 0.00018 |
| 3 | 21.88 | 0.000069 |
Both tests reject decisively, so the residuals are strongly serially correlated: there is systematic, predictable structure left in the error that the model has not captured.
The implication for part (a) is that the reported standard error of 0.085 is not a valid estimate of the sampling uncertainty. Positive serial correlation means the conventional formula treats 42 highly dependent observations as if they were 42 independent ones, so it understates the true standard error and overstates the \(t\)-statistic. The 10.23 in part (a) is not evidence of anything until the inference is repaired.
| Standard error | Estimate | Std. error | \(t\) |
|---|---|---|---|
| Classical OLS | 0.872 | 0.085 | 10.23 |
| Newey–West (3 lags) | 0.872 | 0.118 | 7.40 |
The co-author’s suggestion does what it promises. The estimate is untouched, the standard error rises by about 38%, and the honest \(t\)-statistic is 7.40 rather than 10.23. The conclusion in (a) therefore does not change: population still looks overwhelmingly significant.
That is exactly why HAC is the wrong repair here. It fixes the inference for the equation you wrote down, and nothing else. It does not tell you whether that equation is the right one, and it does not remove the serial correlation, it only prices it. Since the Breusch–Godfrey test is telling us that lags are missing, what we want is a model that contains them, which then delivers both valid inference and the dynamic multipliers that a static regression cannot express. HAC is the right tool when you believe the specification and only distrust the error structure.
Adding inv\(_{t-1}\) and pop\(_{t-1}\) improves matters but does not finish the job: the Breusch–Godfrey test at order one still rejects, with \(p=0.0017\).
The lesson is that “add a lag” is not a ritual but a diagnostic loop. The test tells you whether the dynamics you have written down are rich enough, and here it says they are not. Reporting this specification would leave you in the same position as in part (a), with invalid standard errors, only with more regressors.
With two lags of both variables the residuals are finally clean:
| Term | Estimate | Std. error | \(t\) |
|---|---|---|---|
| Intercept | -33631.7 | 34165.6 | -0.98 |
pop |
4.135 | 9.961 | 0.42 |
pop\(_{t-1}\) |
-10.132 | 18.093 | -0.56 |
pop\(_{t-2}\) |
6.550 | 9.806 | 0.67 |
inv\(_{t-1}\) |
0.864 | 0.148 | 5.85 |
inv\(_{t-2}\) |
-0.520 | 0.152 | -3.42 |
Breusch–Godfrey now gives \(p=0.161\) at order one and \(p=0.351\) at order three, so the model is dynamically complete and its standard errors can be read as usual.
It is tempting to conclude from the table that population does not matter after all, since not one of its three coefficients comes close to significance. That would be the wrong inference. The three population terms are almost perfectly collinear with one another, which inflates each individual standard error while leaving their sum well determined. Testing them jointly gives
\[ F = 5.80, \qquad p = 0.0026, \]
so population does belong in the equation. The long-run multiplier is
\[ \hat\theta = \frac{4.135 - 10.132 + 6.550}{1 - 0.864 + 0.520} = \frac{0.552}{0.656} = 0.842 . \]
Reconciling this with part (a) is the point of the question. The static regression reported 0.872 and the dynamically complete model reports a long-run multiplier of 0.842, which is almost the same number. The static regression was not wrong about the long-run magnitude; it was wrong about everything else. It attached an invalid standard error to that number, it could not distinguish the immediate effect from the eventual one, and it gave no way of knowing whether the figure was a long-run relationship or an artefact. Only the second model licenses the claim.
ardl22 <- lm(inv ~ pop + pop_1 + pop_2 + inv_1 + inv_2, data = housing)
summary(ardl22)
bgtest(ardl22, order = 1)
bgtest(ardl22, order = 3)
# Joint test that the three population terms are all zero
car::linearHypothesis(ardl22, c("pop = 0", "pop_1 = 0", "pop_2 = 0"))
# Long-run multiplier
b <- coef(ardl22)
sum(b[c("pop", "pop_1", "pop_2")]) / (1 - sum(b[c("inv_1", "inv_2")]))ardl22 = smf.ols(
"inv ~ pop + pop_1 + pop_2 + inv_1 + inv_2", data=housing
).fit()
print(ardl22.summary())
print(acorr_breusch_godfrey(ardl22, nlags=1)[:2])
print(acorr_breusch_godfrey(ardl22, nlags=3)[:2])
# Joint test that the three population terms are all zero
print(ardl22.f_test("pop = 0, pop_1 = 0, pop_2 = 0"))
# Long-run multiplier
b = ardl22.params
print((b["pop"] + b["pop_1"] + b["pop_2"]) / (1 - b["inv_1"] - b["inv_2"]))No. The referee has conflated two different failures.
Under strict exogeneity, \(E(u_t\mid X_1,\dots,X_T)=0\), OLS is unbiased in finite samples even when the errors are serially correlated. Serial correlation is a statement about \(\operatorname{Var}(u)\), and unbiasedness does not depend on the shape of the error variance matrix at all.
Under the weaker contemporaneous exogeneity, \(E(u_t\mid X_t)=0\), together with weak dependence, OLS is still consistent, because the sample moments it relies on still converge to their population counterparts.
So a rejected Breusch–Godfrey test on a static regression is not grounds for discarding the estimates. It is grounds for discarding the standard errors, and for asking whether the specification is complete.
Inference. The conventional variance formula \(\hat\sigma^2(X'X)^{-1}\) is derived under the assumption that the errors are uncorrelated across observations, and that assumption has just failed.
With positive serial correlation, which is the usual case in economic time series, the formula treats dependent observations as if each contributed fresh information. It therefore understates the true standard error. Since the \(t\)-statistic divides by that standard error, it is correspondingly overstated, and so tests reject far too often: a nominal 5% test can have a true size several times larger. Confidence intervals are too narrow by the same token. Part (c) of Question 4 is a concrete example: the honest standard error was 38% larger than the reported one.
Take \(Y_t=\alpha+\rho Y_{t-1}+u_t\) with \(u_t=\lambda u_{t-1}+v_t\), \(v_t\) white noise and \(\lambda\neq 0\).
The regressor is \(Y_{t-1}\), and by the model’s own equation \(Y_{t-1}=\alpha+\rho Y_{t-2}+u_{t-1}\), so \(Y_{t-1}\) contains \(u_{t-1}\). The error is \(u_t=\lambda u_{t-1}+v_t\), so \(u_t\) also contains \(u_{t-1}\). Hence
\[ \operatorname{Cov}(Y_{t-1},u_t) = \operatorname{Cov}\big(Y_{t-1},\ \lambda u_{t-1}+v_t\big) = \lambda\operatorname{Cov}(Y_{t-1},u_{t-1}), \]
because \(v_t\) is white noise and therefore uncorrelated with everything dated \(t-1\) or earlier. The remaining covariance is not zero: substituting \(Y_{t-1}=\alpha+\rho Y_{t-2}+u_{t-1}\) gives
\[ \operatorname{Cov}(Y_{t-1},u_{t-1}) = \operatorname{Var}(u_{t-1}) + \rho\operatorname{Cov}(Y_{t-2},u_{t-1}), \]
whose first term is strictly positive. So for any \(\lambda\neq 0\) we have \(\operatorname{Cov}(Y_{t-1},u_t)\neq 0\): the regressor is correlated with the error.
Now even contemporaneous exogeneity fails, so OLS is not merely inefficient but inconsistent: \(\hat\rho\) converges to the wrong number, biased towards whatever value lets the lagged dependent variable absorb the persistence in the error. This is why serial correlation in a dynamic model is a specification problem rather than a standard-error problem, and why HAC standard errors would not rescue it: they would give you an honest confidence interval around an inconsistent estimate.
No, and the reason is precisely what the diagnostics showed. The argument in (c) requires two ingredients: a lagged dependent variable and serial correlation remaining in the error. The specification in Question 4(e) has the first but not the second. Its Breusch–Godfrey tests do not reject, at \(p=0.161\) and \(p=0.351\), so there is no error persistence left for the lagged inv terms to be correlated with.
For those estimates to be untrustworthy in the sense of (c), we would need evidence of remaining autocorrelation, and we tested for it and found none. This is the practical payoff of chasing serial correlation out of the model rather than insuring against it: once the residuals are white noise, the lagged dependent variable is no longer a threat to consistency.
In a long-run equilibrium nothing changes, so \(Y_t=Y_{t-1}=Y_{t-2}=Y_{eq}\) and \(X_t=X_{t-1}=X_{eq}\), and \(E(u_t)=0\). Substituting into the model,
\[ Y_{eq} = 10 + 0.5Y_{eq} + 0.2Y_{eq} + 2.0X_{eq} - 0.8X_{eq}. \]
Collect the \(Y_{eq}\) terms on the left:
\[ Y_{eq}(1 - 0.5 - 0.2) = 10 + (2.0 - 0.8)X_{eq}, \]
so that
\[ Y_{eq} = \frac{10}{1-\sum\phi_i} + \underbrace{\frac{\sum_j\beta_j}{1-\sum_i\phi_i}}_{\theta}X_{eq}. \]
The long-run multiplier is the coefficient on \(X_{eq}\), \(\theta=\sum_j\beta_j/(1-\sum_i\phi_i)\).
\[ \sum_j\beta_j = 2.0 - 0.8 = 1.2, \qquad \sum_i\phi_i = 0.5 + 0.2 = 0.7, \]
\[ \theta = \frac{1.2}{1-0.7} = \frac{1.2}{0.3} = 4.0 . \]
A permanent one-unit increase in \(X\) raises \(Y\) by 4.0 units once all dynamic adjustment is complete.
The contrast with the impact multiplier is the substance of the model. On impact the effect is \(\beta_0=2.0\), exactly half of the eventual response. Iterating the model forward under a permanent unit rise in \(X\) traces the rest of the path:
\[ 2.00,\quad 2.20,\quad 2.70,\quad 2.99,\quad \dots \longrightarrow 4.00 . \]
For instance the second value is \(0.5(2.00)+2.0-0.8=2.20\) and the third is \(0.5(2.20)+0.2(2.00)+1.2=2.70\). The approach is monotone: the negative coefficient on \(X_{t-1}\) slows the buildup, because only \(\beta_0+\beta_1=1.2\) of direct push arrives in each period after the first, but it does not reverse it. A reversal would require a temporary change in \(X\), not a permanent one.
Reporting “the effect of \(X\) on \(Y\) is 2” and “the effect of \(X\) on \(Y\) is 4” are therefore both defensible, provided you say over what horizon.
The numerator is only the direct push that \(X\) gives to \(Y\). But \(Y\) feeds back on itself: whatever change occurs this period is carried into next period with weight \(\sum\phi_i\), and that carry-over is itself carried forward, and so on. The total is a geometric series,
\[ \sum_j\beta_j\left(1 + \textstyle\sum_i\phi_i + (\sum_i\phi_i)^2 + \cdots\right) = \frac{\sum_j\beta_j}{1-\sum_i\phi_i}, \]
which is where the denominator comes from. Equivalently, \(1-\sum\phi_i\) is the fraction of \(Y\) that is genuinely new each period, and dividing by it scales the direct effect up into the accumulated one.
As \(\sum\phi_i\to 1\) the denominator goes to zero and \(\theta\) diverges. Two things should make you cautious. First, the estimate becomes extremely sensitive to \(\sum\hat\phi_i\): with \(\sum\hat\phi_i=0.95\) the denominator is 0.05, so a change in the third decimal place moves \(\theta\) substantially. Second, \(\sum\phi_i=1\) is a unit root in \(Y\), and then there is no long run to speak of, so the quantity you are estimating does not exist. A large long-run multiplier driven by a persistence estimate near one is not a strong finding; it is a warning.
Start from the model and take expectations of both sides:
\[ E(Y_t) = E(\alpha + \rho Y_{t-1} + u_t) = \alpha + \rho E(Y_{t-1}) + E(u_t) = \alpha + \rho E(Y_{t-1}), \]
using \(E(\alpha)=\alpha\) and \(E(u_t)=0\).
Stationarity means the mean does not depend on \(t\), so \(E(Y_t)=E(Y_{t-1})=\mu\). Substituting,
\[ \mu = \alpha + \rho\mu \quad\Longrightarrow\quad \mu(1-\rho) = \alpha \quad\Longrightarrow\quad \mu = \frac{\alpha}{1-\rho}. \]
What breaks at \(\rho=1\). The final step divides by \(1-\rho\), which is zero, so \(\mu\) is not defined. The failure is not merely algebraic. The step before it assumed \(E(Y_t)=E(Y_{t-1})\), and a random walk does not satisfy that: with drift its mean grows without bound, and even without drift its variance \(t\sigma_u^2\) does, so the process has no stationary distribution whose mean we could be computing. The formula does not have a removable singularity at \(\rho=1\); the object it describes ceases to exist.
Substituting the definition of the random walk into the difference gives
\[ \Delta Y_t = (Y_{t-1} + u_t) - Y_{t-1} = u_t, \]
so the first difference is the white noise process. The three conditions follow immediately.
Constant mean. \(E(\Delta Y_t) = E(u_t) = 0\) for every \(t\).
Constant variance. \(\operatorname{Var}(\Delta Y_t) = \operatorname{Var}(u_t) = \sigma^2_u\) for every \(t\). This is the condition the level violated, since \(\operatorname{Var}(Y_t)=t\sigma^2_u\) grows without bound.
Autocovariance depending only on the lag. For \(k>0\), \(\operatorname{Cov}(\Delta Y_t,\Delta Y_{t-k}) = \operatorname{Cov}(u_t,u_{t-k}) = 0\) by the definition of white noise, which is a function of \(k\) alone and not of \(t\).
All three hold, so \(\Delta Y_t\) is covariance stationary.
Why the same remedy is wrong for a trend-stationary series. Suppose instead \(Y_t=\alpha+\delta t+u_t\) with \(u_t\) stationary. Here the trend is deterministic, and including \(t\) as a regressor removes it and leaves the stationary \(u_t\) behind. Differencing gives
\[ \Delta Y_t = \delta + u_t - u_{t-1}, \]
which is stationary but now carries a manufactured MA(1) error: the differenced error has a non-zero first autocovariance, \(-\operatorname{Var}(u)\) in the white noise case, that was not there before. We have introduced serial correlation and thrown away information for no gain. The general rule is that a stochastic trend needs differencing and a deterministic trend needs a trend term, which is why the choice is a testable question rather than a matter of habit.
The loop below re-estimates the AR(2) model once for every year of the evaluation period, always on data up to \(t-1\) only, and stores the resulting one-step-ahead forecast alongside the two benchmarks. The important discipline is that nothing inside the loop is allowed to see year \(t\): the training sample, the estimated coefficients and the training mean all stop at \(t-1\).
library(dplyr)
oos_years <- (max(fertil$year) - 19):max(fertil$year)
oos <- lapply(oos_years, function(yr) {
train <- fertil |> filter(year < yr)
test <- fertil |> filter(year == yr)
ar2 <- lm(gfr ~ gfr_1 + gfr_2, data = train)
data.frame(
year = yr,
actual = test$gfr,
ar2 = unname(predict(ar2, newdata = test)),
random_walk = test$gfr_1,
train_mean = mean(train$gfr)
)
}) |>
bind_rows()
head(oos)fertil["gfr_2"] = fertil["gfr"].shift(2)
oos_years = range(int(fertil["year"].max()) - 19,
int(fertil["year"].max()) + 1)
rows = []
for yr in oos_years:
train = fertil[fertil["year"] < yr]
test = fertil[fertil["year"] == yr]
ar2 = smf.ols("gfr ~ gfr_1 + gfr_2", data=train).fit()
rows.append({
"year": yr,
"actual": test["gfr"].iloc[0],
"ar2": ar2.predict(test).iloc[0],
"random_walk": test["gfr_1"].iloc[0],
"train_mean": train["gfr"].mean(),
})
oos = pd.DataFrame(rows)
print(oos.head())| Rule | RMSE | MAE |
|---|---|---|
| AR(2) | 3.156 | 2.200 |
| Random walk, \(\hat Y_t = Y_{t-1}\) | 3.611 | 2.590 |
| Training-sample mean | 26.933 | 25.597 |
The AR(2) forecasts best, and here the ranking is the same under both metrics: it beats the random walk by about 13% on RMSE and 15% on MAE. That agreement is worth noticing rather than assuming, because RMSE squares the errors and MAE does not, so a handful of large misses can reverse the ordering. When two metrics disagree, you have to say which loss function you are minimising before you declare a winner.
The margin over the random walk is also modest. A series this persistent is largely predicted by its own last value, so most of what any model can achieve is already achieved by the naive benchmark. Reporting the benchmark alongside the model is what makes the improvement interpretable.
Because the series is not mean-reverting on any relevant timescale. Fertility falls almost monotonically from the baby boom peak through the whole forecast period, so the mean of the training sample always sits far above the current level, and it lags further behind as the decline continues. By 1984 the training mean is about 96 while the actual value is 65.4, an error of over 30.
This is a warning about long-horizon forecasts from a stationary model. In Question 1(e) we saw the AR(2) forecast drifting towards its unconditional mean of 77.5, and the sample-mean rule is simply the limiting case of that behaviour, applied immediately. If the unconditional mean is not a credible statement about where the series is going, then neither is any forecast that converges to it. The short-horizon forecasts in the table above are good precisely because they lean on the recent past rather than on the estimated mean, and their quality should not be extrapolated to horizons where the mean takes over.