import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as smDataset One: GDP¶
The following data is on US GDP from FRED (https://
gdp = pd.read_csv("GDP_29Oct2026.csv")
print(gdp.head())
y = gdp['GDP'].to_numpy()
plt.plot(y, color = 'black')
plt.xlabel('Quarter')
plt.ylabel('Billions of Dollars')
plt.title('Gross Domestic Product')
plt.show() observation_date GDP
0 1947-01-01 243.164
1 1947-04-01 245.968
2 1947-07-01 249.585
3 1947-10-01 259.745
4 1948-01-01 265.742
We will fit AR(p) models to this dataset, and use the fitted models to predict future values. Let us leave the last 50 observations (corresponding to 12.5 years; note that this is quarterly data) as test values, and the remaining data form the training data. We will fit AR(p) models to the training data, predict the test values and then compare the predictions with the actual test values.
n = len(y)
tme = range(1, n+1)
n_test = 50
#n_test = 12
n_train = n - n_test
y_train = y[:n_train]
tme_train = tme[:n_train]
y_test = y[n_train:]
tme_test = tme[n_train:]plt.plot(tme_train, y_train, color='black')
plt.xlabel('Quarter')
plt.ylabel('Billions of Dollars')
plt.title('Gross Domestic Product — Training Data')
plt.show()
Below, we fit AR(1) to the training data. There are two ways of doing this. Either manually create and and use OLS, or use the AutoReg function from statsmodels.
p = 1
yreg = y_train[p:] #these are the response values in the autoregression
Xmat = np.ones((n_train-p, 1)) #this will be the design matrix (X) in the autoregression
for j in range(1, p+1):
col = y_train[p-j : n_train-j].reshape(-1, 1)
Xmat = np.column_stack([Xmat, col])
armod = sm.OLS(yreg, Xmat).fit()
print(armod.params)
print(armod.summary())[19.18439505 1.0087699 ]
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 1.000
Model: OLS Adj. R-squared: 1.000
Method: Least Squares F-statistic: 2.455e+06
Date: Fri, 30 Oct 2026 Prob (F-statistic): 0.00
Time: 21:36:28 Log-Likelihood: -1447.8
No. Observations: 268 AIC: 2900.
Df Residuals: 266 BIC: 2907.
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 19.1844 4.616 4.156 0.000 10.096 28.272
x1 1.0088 0.001 1566.951 0.000 1.008 1.010
==============================================================================
Omnibus: 215.657 Durbin-Watson: 1.002
Prob(Omnibus): 0.000 Jarque-Bera (JB): 5540.618
Skew: -2.941 Prob(JB): 0.00
Kurtosis: 24.484 Cond. No. 1.01e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.01e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
from statsmodels.tsa.ar_model import AutoReg
armod_sm = AutoReg(y_train, lags = p).fit()
print(armod_sm.summary()) AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 269
Model: AutoReg(1) Log Likelihood -1447.753
Method: Conditional MLE S.D. of innovations 53.685
Date: Fri, 30 Oct 2026 AIC 2901.507
Time: 21:36:28 BIC 2912.279
Sample: 1 HQIC 2905.833
269
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 19.1844 4.598 4.172 0.000 10.172 28.197
y.L1 1.0088 0.001 1572.831 0.000 1.008 1.010
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9913 +0.0000j 0.9913 0.0000
-----------------------------------------------------------------------------
Here are some observations comparing the above two outputs:
The parameter estimates are exactly the same.
The regression summary gives t-scores corresponding to each parameter estimate while the AutoReg summary only gives z-scores
The standard errors are slightly different.
The standard errors given by the regression summary correspond to the square roots of the diagonal entries of:
while the standard errors reported by AutoReg summary correspond to the square roots of the diagonal entries of:
Here is how the fitted model can be used to predict the test observations.
k = n_test
fcast = armod_sm.get_prediction(start = n_train, end = n_train+k-1)
fcast_mean = fcast.predicted_mean #this gives point predictions for the future values of the time series The method that the above function uses to calculate the predictions was detailed in class; and is captured in the code below.
#Predictions
yhat = np.concatenate([y_train.astype(float), np.full(k, -9999)]) #extend data by k placeholder values
phi_vals = armod_sm.params
for i in range(1, k+1):
ans = phi_vals[0]
for j in range(1, p+1):
ans += phi_vals[j] * yhat[n_train+i-j-1]
yhat[n_train+i-1] = ans
predvalues = yhat[n_train:]
#Check that both predictions are identical:
print(np.column_stack([predvalues, fcast_mean]))
preds_nolog = fcast_mean[[17367.744854 17367.744854 ]
[17539.24265126 17539.24265126]
[17712.24446723 17712.24446723]
[17886.76349198 17886.76349198]
[18062.8130313 18062.8130313 ]
[18240.40650764 18240.40650764]
[18419.55746118 18419.55746118]
[18600.27955083 18600.27955083]
[18782.58655532 18782.58655532]
[18966.49237418 18966.49237418]
[19152.01102885 19152.01102885]
[19339.15666375 19339.15666375]
[19527.94354734 19527.94354734]
[19718.38607319 19718.38607319]
[19910.49876114 19910.49876114]
[20104.29625833 20104.29625833]
[20299.79334038 20299.79334038]
[20497.00491247 20497.00491247]
[20695.94601052 20695.94601052]
[20896.6318023 20896.6318023 ]
[21099.07758859 21099.07758859]
[21303.29880438 21303.29880438]
[21509.31102 21509.31102 ]
[21717.12994235 21717.12994235]
[21926.77141607 21926.77141607]
[22138.25142475 22138.25142475]
[22351.58609216 22351.58609216]
[22566.79168347 22566.79168347]
[22783.8846065 22783.8846065 ]
[23002.88141297 23002.88141297]
[23223.79879974 23223.79879974]
[23446.65361011 23446.65361011]
[23671.4628351 23671.4628351 ]
[23898.24361472 23898.24361472]
[24127.01323932 24127.01323932]
[24357.78915087 24357.78915087]
[24590.58894431 24590.58894431]
[24825.43036888 24825.43036888]
[25062.33132948 25062.33132948]
[25301.30988805 25301.30988805]
[25542.3842649 25542.3842649 ]
[25785.57284016 25785.57284016]
[26030.89415515 26030.89415515]
[26278.36691376 26278.36691376]
[26528.00998396 26528.00998396]
[26779.84239916 26779.84239916]
[27033.8833597 27033.8833597 ]
[27290.1522343 27290.1522343 ]
[27548.66856155 27548.66856155]
[27809.45205139 27809.45205139]]
Below we plot the predictions along with the actual test observations.
plt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, fcast_mean, label = 'Forecast', color = 'black')
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
Next we discuss prediction uncertainty quantification. The prediction standard errors, and the associated uncertainty intervals, are calculated as follows.
fcast_se = fcast.se_mean
alpha = 0.05
from scipy import stats
z_alpha_half = stats.norm.ppf(1 - alpha/2)
predlower = predvalues - z_alpha_half * fcast_se
predupper = predvalues + z_alpha_half * fcast_se
fcast_int = fcast.conf_int()
print(np.column_stack([predlower, predupper, fcast_int]))
[[17262.525004 17472.96470399 17262.525004 17472.96470399]
[17389.78539276 17688.69990976 17389.78539276 17688.69990976]
[17528.39071614 17896.09821831 17528.39071614 17896.09821831]
[17673.52910847 18099.99787549 17673.52910847 18099.99787549]
[17823.35294746 18302.27311514 17823.35294746 18302.27311514]
[17976.92498864 18503.88802663 17976.92498864 18503.88802663]
[18133.69608893 18705.41883342 18133.69608893 18705.41883342]
[18293.31453482 18907.24456685 18293.31453482 18907.24456685]
[18455.54111949 19109.63199114 18455.54111949 19109.63199114]
[18620.20608372 19312.77866463 18620.20608372 19312.77866463]
[18787.1851721 19516.83688561 18787.1851721 19516.83688561]
[18956.385361 19721.9279665 18956.385361 19721.9279665 ]
[19127.73587643 19928.15121825 19127.73587643 19928.15121825]
[19301.18228807 20135.58985832 19301.18228807 20135.58985832]
[19476.68248368 20344.3150386 19476.68248368 20344.3150386 ]
[19654.20384085 20554.3886758 19654.20384085 20554.3886758 ]
[19833.72118831 20765.86549245 19833.72118831 20765.86549245]
[20015.21530335 20978.79452159 20015.21530335 20978.79452159]
[20198.67178301 21193.22023803 20198.67178301 21193.22023803]
[20384.0801813 21409.1834233 20384.0801813 21409.1834233 ]
[20571.43334007 21626.72183711 20571.43334007 21626.72183711]
[20760.72686292 21845.87074584 20760.72686292 21845.87074584]
[20951.95869656 22066.66334344 20951.95869656 22066.66334344]
[21145.12879424 22289.13109046 21145.12879424 22289.13109046]
[21340.23884233 22513.30398981 21340.23884233 22513.30398981]
[21537.29203624 22739.21081325 21537.29203624 22739.21081325]
[21736.29289545 22966.87928886 21736.29289545 22966.87928886]
[21937.24710942 23196.33625752 21937.24710942 23196.33625752]
[22140.16140861 23427.6078044 22140.16140861 23427.6078044 ]
[22345.04345566 23660.71937028 22345.04345566 23660.71937028]
[22551.90175316 23895.69584632 22551.90175316 23895.69584632]
[22760.74556499 24132.56165524 22760.74556499 24132.56165524]
[22971.58484895 24371.34082125 22971.58484895 24371.34082125]
[23184.43019882 24612.05703063 23184.43019882 24612.05703063]
[23399.29279424 24854.73368441 23399.29279424 24854.73368441]
[23616.18435726 25099.39394448 23616.18435726 25099.39394448]
[23835.11711456 25346.06077406 23835.11711456 25346.06077406]
[24056.10376435 25594.75697341 24056.10376435 25594.75697341]
[24279.15744752 25845.50521145 24279.15744752 25845.50521145]
[24504.29172214 26098.32805395 24504.29172214 26098.32805395]
[24731.52054119 26353.24798862 24731.52054119 26353.24798862]
[24960.85823272 26610.28744761 24960.85823272 26610.28744761]
[25192.31948246 26869.46882783 25192.31948246 26869.46882783]
[25425.91931834 27130.81450919 25425.91931834 27130.81450919]
[25661.67309677 27394.34687116 25661.67309677 27394.34687116]
[25899.59649047 27660.08830785 25899.59649047 27660.08830785]
[26139.70547763 27928.06124177 26139.70547763 27928.06124177]
[26382.0163323 28198.28813631 26382.0163323 28198.28813631]
[26626.54561575 28470.79150735 26626.54561575 28470.79150735]
[26873.31016892 28745.59393386 26873.31016892 28745.59393386]]
Below we plot the predictions along with uncertainty intervals.
#Plotting predictions along with uncertainty:
plt.plot(tme_train, y_train, label = 'Original Data')
plt.plot(tme_test, predvalues, label = 'Forecast', color = 'black')
plt.plot(tme_test, predlower, color = 'green', label = 'Prediction lower bound')
plt.plot(tme_test, predupper, color = 'green', label = 'Prediction upper bound')
plt.plot(tme_test, y_test, color = 'red', label = 'Actual Future Values')
plt.legend()
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.show()
In this example, the predicted values are somewhat below the actual test values.
Instead of applying AR models, directly to the raw data, it is common practice to apply them to the logarithms.
For GDP (and also for the house prices data which we shall look at next), it is often natural to think of fluctuations in percentage terms. A 5% fluctuation is twice as large in dollars when the level of the series doubles, but its size on the log scale is unchanged: multiplying by 1.05 adds to . Taking logarithms can therefore make the residual spread more stable, better matching the constant error variance assumption of an AR model.
ylog_train = np.log(y_train)
ylog_test = np.log(y_test)armod_sm = AutoReg(ylog_train, lags = p).fit()
print(armod_sm.summary()) AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 269
Model: AutoReg(1) Log Likelihood 830.433
Method: Conditional MLE S.D. of innovations 0.011
Date: Fri, 30 Oct 2026 AIC -1654.866
Time: 21:40:23 BIC -1644.093
Sample: 1 HQIC -1650.539
269
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0283 0.004 7.232 0.000 0.021 0.036
y.L1 0.9984 0.000 2014.926 0.000 0.997 0.999
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.0016 +0.0000j 1.0016 0.0000
-----------------------------------------------------------------------------
Predictions are obtained as follows.
fcast = armod_sm.get_prediction(start = n_train, end = n_train+k-1)
fcast_mean = fcast.predicted_mean #this gives the point predictionsPrediction uncertainty intervals are obtained as follows.
fcast_se = fcast.se_mean
alpha = 0.05
z_alpha_half = stats.norm.ppf(1 - alpha/2)
predlower = fcast_mean - z_alpha_half * fcast_se
predupper = fcast_mean + z_alpha_half * fcast_se
fcast_int = fcast.conf_int()
print(np.column_stack([predlower, predupper, fcast_int]))[[ 9.74390031 9.78668836 9.74390031 9.78668836]
[ 9.74780358 9.80826671 9.74780358 9.80826671]
[ 9.7537592 9.82775201 9.7537592 9.82775201]
[ 9.76077004 9.84614143 9.76077004 9.84614143]
[ 9.76844957 9.8638216 9.76844957 9.8638216 ]
[ 9.77659938 9.88099097 9.77659938 9.88099097]
[ 9.78510153 9.89776754 9.78510153 9.89776754]
[ 9.79387913 9.9142283 9.79387913 9.9142283 ]
[ 9.80287871 9.93042675 9.80287871 9.93042675]
[ 9.81206134 9.94640189 9.81206134 9.94640189]
[ 9.82139761 9.96218321 9.82139761 9.96218321]
[ 9.83086465 9.97779363 9.83086465 9.97779363]
[ 9.84044427 9.99325141 9.84044427 9.99325141]
[ 9.85012171 10.00857138 9.85012171 10.00857138]
[ 9.85988481 10.02376575 9.85988481 10.02376575]
[ 9.86972339 10.03884477 9.86972339 10.03884477]
[ 9.87962885 10.0538171 9.87962885 10.0538171 ]
[ 9.88959383 10.06869017 9.88959383 10.06869017]
[ 9.89961198 10.0834704 9.89961198 10.0834704 ]
[ 9.90967777 10.09816336 9.90967777 10.09816336]
[ 9.91978637 10.11277397 9.91978637 10.11277397]
[ 9.92993349 10.12730657 9.92993349 10.12730657]
[ 9.94011534 10.141765 9.94011534 10.141765 ]
[ 9.95032855 10.15615272 9.95032855 10.15615272]
[ 9.96057008 10.17047282 9.96057008 10.17047282]
[ 9.97083719 10.1847281 9.97083719 10.1847281 ]
[ 9.98112741 10.19892109 9.98112741 10.19892109]
[ 9.9914385 10.2130541 9.9914385 10.2130541 ]
[10.00176842 10.22712924 10.00176842 10.22712924]
[10.01211529 10.24114843 10.01211529 10.24114843]
[10.02247741 10.25511346 10.02247741 10.25511346]
[10.03285319 10.26902595 10.03285319 10.26902595]
[10.04324119 10.28288743 10.04324119 10.28288743]
[10.05364006 10.2966993 10.05364006 10.2966993 ]
[10.06404856 10.31046286 10.06404856 10.31046286]
[10.07446554 10.32417933 10.07446554 10.32417933]
[10.08488992 10.33784985 10.08488992 10.33784985]
[10.0953207 10.35147547 10.0953207 10.35147547]
[10.10575696 10.36505718 10.10575696 10.36505718]
[10.11619781 10.37859593 10.11619781 10.37859593]
[10.12664245 10.39209259 10.12664245 10.39209259]
[10.13709009 10.405548 10.13709009 10.405548 ]
[10.14754003 10.41896292 10.14754003 10.41896292]
[10.15799159 10.43233811 10.15799159 10.43233811]
[10.16844412 10.44567425 10.16844412 10.44567425]
[10.17889703 10.45897202 10.17889703 10.45897202]
[10.18934974 10.47223204 10.18934974 10.47223204]
[10.19980172 10.48545491 10.19980172 10.48545491]
[10.21025246 10.4986412 10.21025246 10.4986412 ]
[10.22070148 10.51179145 10.22070148 10.51179145]]
plt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, np.exp(fcast_mean), label = 'Forecast', color = 'black') #Note the exponentiation
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.plot(tme_test, np.exp(predlower), color = 'green', label = 'Prediction lower bound')
plt.plot(tme_test, np.exp(predupper), color = 'green', label = 'Prediction upper bound')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
Note that now the predictions are slightly closer to the actual observations.
plt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, np.exp(fcast_mean), label = 'Predictions with logs', color = 'black') #Note the exponentiation
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.plot(tme_test, preds_nolog, color = 'darkgreen', label = 'Predictions without logs')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
Residuals: raw data versus logarithms¶
Below we plot the residuals from the AR fit for both the models (AR fit to the raw data and AR fit to the logarithms).
resid_raw = AutoReg(y_train, lags=p).fit().resid
resid_log = armod_sm.resid
fig, axes = plt.subplots(2, 1, figsize=(11, 6), sharex=True)
for ax, resid, label in zip(axes, [resid_raw, resid_log],
['Raw data (billions of dollars)', 'Log data (log units)']):
ax.plot(tme_train[p:], resid, color='black', linewidth=1)
ax.axhline(0, color='red', linestyle='--')
ax.set(title=label, ylabel='Residual')
axes[-1].set_xlabel('Quarter')
plt.tight_layout()
plt.show()
The raw residuals become much more variable as spending grows. Taking logarithms substantially reduces this pattern, making the constant error variance assumption more plausible. Compare how the spread changes over time, not the numerical sizes of residuals on these different scales.
Dataset Two: House Price Data from FRED¶
Consider the following quarterly dataset (again from FRED: https://
hprice = pd.read_csv('ASPUS_October2026.csv')
print(hprice.head())
y = hprice['ASPUS'].to_numpy()
plt.figure(figsize = (12, 6))
plt.plot(y)
plt.xlabel('Year')
plt.ylabel('Dollars')
plt.title('Average Sales Price of Houses Sold in the US')
plt.show() observation_date ASPUS
0 1963-01-01 19300
1 1963-04-01 19400
2 1963-07-01 19200
3 1963-10-01 19600
4 1964-01-01 19600
Training-test split:
n = len(y)
tme = range(1, n+1)
#n_test = 50
n_test = 100
n_train = n - n_test
y_train = y[:n_train]
tme_train = tme[:n_train]
y_test = y[n_train:]
tme_test = tme[n_train:]Below we fit AR(1) to the raw data.
p = 1
from statsmodels.tsa.ar_model import AutoReg
armod_sm = AutoReg(y_train, lags = p).fit()
print(armod_sm.summary()) AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 155
Model: AutoReg(1) Log Likelihood -1421.301
Method: Conditional MLE S.D. of innovations 2465.844
Date: Fri, 30 Oct 2026 AIC 2848.602
Time: 21:55:35 BIC 2857.713
Sample: 1 HQIC 2852.303
155
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 671.9911 368.218 1.825 0.068 -49.702 1393.684
y.L1 1.0059 0.003 302.244 0.000 0.999 1.012
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9941 +0.0000j 0.9941 0.0000
-----------------------------------------------------------------------------
Below we obtain predictions.
k = n_test
fcast = armod_sm.get_prediction(start = n_train, end = n_train+k-1)
fcast_mean = fcast.predicted_mean #this gives point predictions for the future values of the time series
preds_nolog = fcast_meanplt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, fcast_mean, label = 'Forecast (no logs)', color = 'black')
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
fcast_se = fcast.se_mean
alpha = 0.05
z_alpha_half = stats.norm.ppf(1 - alpha/2)
predlower = fcast_mean - z_alpha_half * fcast_se
predupper = fcast_mean + z_alpha_half * fcast_se
fcast_int = fcast.conf_int()
print(np.column_stack([predlower, predupper, fcast_int]))[[204870.58737551 214536.51969315 204870.58737551 214536.51969315]
[204763.25865112 218473.51896037 204763.25865112 218473.51896037]
[205123.821886 221965.32346879 205123.821886 221965.32346879]
[205729.75702797 225234.58778946 205729.75702797 225234.58778946]
[206495.17430941 228367.3370061 206495.17430941 228367.3370061 ]
[207376.18106667 231407.59990046 207376.18106667 231407.59990046]
[208346.80078729 234381.48990978 208346.80078729 234381.48990978]
[209390.17604703 237306.00219464 209390.17604703 237306.00219464]
[210494.64798737 240192.93416618 210494.64798737 240192.93416618]
[211651.76725522 243050.87455129 211651.76725522 243050.87455129]
[212855.18728792 245886.31011244 212855.18728792 245886.31011244]
[214100.00424011 248704.28572573 214100.00424011 248704.28572573]
[215382.34130984 251508.82005967 215382.34130984 251508.82005967]
[216699.07525994 254303.17905882 216699.07525994 254303.17905882]
[218047.64988628 257090.06248045 218047.64988628 257090.06248045]
[219425.94488395 259871.73503341 219425.94488395 259871.73503341]
[220832.18125492 262650.12097553 220832.18125492 262650.12097553]
[222264.85154358 265426.87388308 222264.85154358 265426.87388308]
[223722.66737955 268203.42911315 223722.66737955 268203.42911315]
[225204.51935796 270981.04392843 225204.51935796 270981.04392843]
[226709.44589015 273760.82865175 226709.44589015 273760.82865175]
[228236.60869166 276543.7711832 228236.60869166 276543.7711832 ]
[229785.27325923 279330.75652839 229785.27325923 279330.75652839]
[231354.79315084 282122.58252368 231354.79315084 282122.58252368]
[232944.5972022 284919.97262499 232944.5972022 284919.97262499]
[234554.17903691 287723.58640289 234554.17903691 287723.58640289]
[236183.08838755 290534.02822692 236183.08838755 290534.02822692]
[237830.9238602 293351.85450644 237830.9238602 293351.85450644]
[239497.32686011 296177.57977035 239497.32686011 296177.57977035]
[241181.97645903 299011.68180524 241181.97645903 299011.68180524]
[242884.5850321 301854.60602394 242884.5850321 301854.60602394]
[244604.8945282 304706.76920069 244604.8945282 304706.76920069]
[246342.67326511 307568.56268157 246342.67326511 307568.56268157]
[248097.71316219 310440.35515732 248097.71316219 310440.35515732]
[249869.82733998 313322.49506938 249869.82733998 313322.49506938]
[251658.84802912 316215.31270665 251658.84802912 316215.31270665]
[253464.6247414 319119.12204003 253464.6247414 319119.12204003]
[255287.022664 322034.22233399 255287.022664 322034.22233399]
[257125.92124464 324960.89956707 257125.92124464 324960.89956707]
[258981.21294074 327899.42768857 258981.21294074 327899.42768857]
[260852.80211001 330850.0697337 260852.80211001 330850.0697337 ]
[262740.60402357 333813.07881639 262740.60402357 333813.07881639]
[264644.54398538 336788.69901564 264644.54398538 336788.69901564]
[266564.55654464 339777.16616914 266564.55654464 339777.16616914]
[268500.58478924 342778.70858567 268500.58478924 342778.70858567]
[270452.57971055 345793.54768639 270452.57971055 345793.54768639]
[272420.49963091 348821.89858328 272420.49963091 348821.89858328]
[274404.30968652 351863.97060236 274404.30968652 351863.97060236]
[276403.98135931 354919.96775792 276403.98135931 354919.96775792]
[278419.49205236 357990.08918318 278419.49205236 357990.08918318]
[280450.82470401 361074.52952239 280450.82470401 361074.52952239]
[282497.96743651 364173.47928837 282497.96743651 364173.47928837]
[284560.91323549 367287.12518924 284560.91323549 367287.12518924]
[286639.65965717 370415.65042755 286639.65965717 370415.65042755]
[288734.20856033 373559.23497459 288734.20856033 373559.23497459]
[290844.56586061 376718.05582239 290844.56586061 376718.05582239]
[292970.74130504 379892.28721568 292970.74130504 379892.28721568]
[295112.74826465 383082.10086562 295112.74826465 383082.10086562]
[297270.6035436 386287.66614716 297270.6035436 386287.66614716]
[299444.32720323 389509.1502815 299444.32720323 389509.1502815 ]
[301633.94239963 392746.71850506 301633.94239963 392746.71850506]
[303839.47523356 396000.53422611 303839.47523356 396000.53422611]
[306060.95461156 399270.75917032 306060.95461156 399270.75917032]
[308298.41211733 402557.553516 308298.41211733 402557.553516 ]
[310551.88189247 405861.07602007 310551.88189247 405861.07602007]
[312821.40052575 409181.48413553 312821.40052575 409181.48413553]
[315107.00695026 412518.93412108 315107.00695026 412518.93412108]
[317408.74234776 415873.58114361 317408.74234776 415873.58114361]
[319726.6500596 419245.57937406 319726.6500596 419245.57937406]
[322060.77550379 422635.08207736 322060.77550379 422635.08207736]
[324411.16609765 426042.24169663 324411.16609765 426042.24169663]
[326777.8711856 429467.20993241 326777.8711856 429467.20993241]
[329160.94197177 432910.13781702 329160.94197177 432910.13781702]
[331560.43145701 436371.1757847 331560.43145701 436371.1757847 ]
[333976.39437998 439850.47373754 333976.39437998 439850.47373754]
[336408.887162 443348.18110782 336408.887162 443348.18110782]
[338857.9678555 446864.44691674 338857.9678555 446864.44691674]
[341323.69609559 450399.41983006 341323.69609559 450399.41983006]
[343806.13305478 453953.2482107 343806.13305478 453953.2482107 ]
[346305.34140045 457526.08016854 346305.34140045 457526.08016854]
[348821.38525499 461118.06360765 348821.38525499 461118.06360765]
[351354.33015831 464729.3462712 351354.33015831 464729.3462712 ]
[353904.24303274 468360.075784 353904.24303274 468360.075784 ]
[356471.19214997 472010.39969306 356471.19214997 472010.39969306]
[359055.2471 475680.46550614 359055.2471 475680.46550614]
[361656.47876199 479370.42072853 361656.47876199 479370.42072853]
[364274.95927683 483080.41289809 364274.95927683 483080.41289809]
[366910.76202129 486810.58961873 366910.76202129 486810.58961873]
[369563.96158382 490561.09859239 369563.96158382 490561.09859239]
[372234.63374169 494332.08764962 372234.63374169 494332.08764962]
[374922.85543948 498123.70477891 374922.85543948 498123.70477891]
[377628.70476893 501936.0981547 377628.70476893 501936.0981547 ]
[380352.26094988 505769.41616434 380352.26094988 505769.41616434]
[383093.60431238 509623.80743391 383093.60431238 509623.80743391]
[385852.81627987 513499.42085304 385852.81627987 513499.42085304]
[388629.97935333 517396.40559884 388629.97935333 517396.40559884]
[391425.1770964 521314.91115883 391425.1770964 521314.91115883]
[394238.49412134 525255.08735309 394238.49412134 525255.08735309]
[397070.01607592 529217.08435566 397070.01607592 529217.08435566]
[399919.82963098 533201.05271509 399919.82963098 533201.05271509]]
#Plotting predictions along with uncertainty:
plt.plot(tme_train, y_train, label = 'Original Data')
plt.plot(tme_test, fcast_mean, label = 'Forecast', color = 'black')
plt.plot(tme_test, predlower, color = 'green', label = 'Prediction lower bound')
plt.plot(tme_test, predupper, color = 'green', label = 'Prediction upper bound')
plt.plot(tme_test, y_test, color = 'red', label = 'Actual Future Values')
plt.legend()
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.show()
Now we take logs.
ylog_train = np.log(y_train)
ylog_test = np.log(y_test)armod_sm = AutoReg(ylog_train, lags = p).fit()
print(armod_sm.summary()) AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 155
Model: AutoReg(1) Log Likelihood 361.891
Method: Conditional MLE S.D. of innovations 0.023
Date: Fri, 30 Oct 2026 AIC -717.783
Time: 21:55:53 BIC -708.672
Sample: 1 HQIC -714.082
155
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0573 0.027 2.139 0.032 0.005 0.110
y.L1 0.9963 0.002 416.930 0.000 0.992 1.001
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.0038 +0.0000j 1.0038 0.0000
-----------------------------------------------------------------------------
fcast = armod_sm.get_prediction(start = n_train, end = n_train+k-1)
fcast_mean = fcast.predicted_mean #this gives the point predictionsplt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, np.exp(fcast_mean), label = 'Predictions with logs', color = 'black') #Note the exponentiation
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.plot(tme_test, preds_nolog, color = 'darkgreen', label = 'Predictions without logs')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
fcast_se = fcast.se_mean
alpha = 0.05
z_alpha_half = stats.norm.ppf(1 - alpha/2)
predlower = fcast_mean - z_alpha_half * fcast_se
predupper = fcast_mean + z_alpha_half * fcast_se
fcast_int = fcast.conf_int()
print(np.column_stack([predlower, predupper, fcast_int]))[[12.21053987 12.30100338 12.21053987 12.30100338]
[12.2033213 12.33101686 12.2033213 12.33101686]
[12.20047263 12.35657513 12.20047263 12.35657513]
[12.19987842 12.37979397 12.19987842 12.37979397]
[12.20071763 12.40149472 12.20071763 12.40149472]
[12.20256838 12.42209958 12.20256838 12.42209958]
[12.20517962 12.44185993 12.20517962 12.44185993]
[12.20838747 12.46093995 12.20838747 12.46093995]
[12.21207794 12.47945394 12.21207794 12.47945394]
[12.21616788 12.49748536 12.21616788 12.49748536]
[12.22059445 12.51509739 12.22059445 12.51509739]
[12.22530874 12.53233922 12.22530874 12.53233922]
[12.2302718 12.54925011 12.2302718 12.54925011]
[12.23545204 12.56586197 12.23545204 12.56586197]
[12.24082335 12.58220121 12.24082335 12.58220121]
[12.24636393 12.59828994 12.24636393 12.59828994]
[12.2520553 12.61414694 12.2520553 12.61414694]
[12.25788165 12.62978832 12.25788165 12.62978832]
[12.26382933 12.64522804 12.26382933 12.64522804]
[12.26988646 12.66047827 12.26988646 12.66047827]
[12.27604262 12.67554973 12.27604262 12.67554973]
[12.28228862 12.69045192 12.28228862 12.69045192]
[12.2886163 12.7051933 12.2886163 12.7051933 ]
[12.29501836 12.71978145 12.29501836 12.71978145]
[12.30148827 12.7342232 12.30148827 12.7342232 ]
[12.30802014 12.74852474 12.30802014 12.74852474]
[12.31460866 12.76269167 12.31460866 12.76269167]
[12.32124898 12.77672912 12.32124898 12.77672912]
[12.32793669 12.7906418 12.32793669 12.7906418 ]
[12.33466779 12.80443402 12.33466779 12.80443402]
[12.34143856 12.81810975 12.34143856 12.81810975]
[12.34824563 12.83167268 12.34824563 12.83167268]
[12.35508586 12.84512622 12.35508586 12.84512622]
[12.36195637 12.85847353 12.36195637 12.85847353]
[12.36885448 12.87171759 12.36885448 12.87171759]
[12.37577771 12.88486115 12.37577771 12.88486115]
[12.38272377 12.8979068 12.38272377 12.8979068 ]
[12.38969049 12.91085696 12.38969049 12.91085696]
[12.39667588 12.92371393 12.39667588 12.92371393]
[12.40367806 12.93647986 12.40367806 12.93647986]
[12.41069528 12.94915676 12.41069528 12.94915676]
[12.41772591 12.96174657 12.41772591 12.96174657]
[12.42476839 12.9742511 12.42476839 12.9742511 ]
[12.43182129 12.98667207 12.43182129 12.98667207]
[12.43888323 12.99901112 12.43888323 12.99901112]
[12.44595293 13.0112698 12.44595293 13.0112698 ]
[12.45302919 13.0234496 12.45302919 13.0234496 ]
[12.46011086 13.03555193 12.46011086 13.03555193]
[12.46719686 13.04757814 12.46719686 13.04757814]
[12.47428617 13.05952952 12.47428617 13.05952952]
[12.48137783 13.07140731 12.48137783 13.07140731]
[12.48847091 13.08321268 12.48847091 13.08321268]
[12.49556454 13.09494678 12.49556454 13.09494678]
[12.50265791 13.10661069 12.50265791 13.10661069]
[12.50975023 13.11820546 12.50975023 13.11820546]
[12.51684075 13.12973209 12.51684075 13.12973209]
[12.52392877 13.14119156 12.52392877 13.14119156]
[12.53101361 13.15258479 12.53101361 13.15258479]
[12.53809464 13.1639127 12.53809464 13.1639127 ]
[12.54517124 13.17517615 12.54517124 13.17517615]
[12.55224283 13.18637597 12.55224283 13.18637597]
[12.55930886 13.19751299 12.55930886 13.19751299]
[12.56636878 13.20858798 12.56636878 13.20858798]
[12.57342211 13.21960171 12.57342211 13.21960171]
[12.58046836 13.23055492 12.58046836 13.23055492]
[12.58750706 13.24144831 12.58750706 13.24144831]
[12.59453778 13.25228258 12.59453778 13.25228258]
[12.60156008 13.26305841 12.60156008 13.26305841]
[12.60857358 13.27377645 12.60857358 13.27377645]
[12.61557788 13.28443733 12.61557788 13.28443733]
[12.62257262 13.29504167 12.62257262 13.29504167]
[12.62955744 13.30559006 12.62955744 13.30559006]
[12.636532 13.31608311 12.636532 13.31608311]
[12.64349599 13.32652136 12.64349599 13.32652136]
[12.65044908 13.33690539 12.65044908 13.33690539]
[12.65739098 13.34723572 12.65739098 13.34723572]
[12.66432142 13.3575129 12.66432142 13.3575129 ]
[12.67124011 13.36773742 12.67124011 13.36773742]
[12.67814679 13.3779098 12.67814679 13.3779098 ]
[12.68504121 13.38803053 12.68504121 13.38803053]
[12.69192313 13.39810008 12.69192313 13.39810008]
[12.69879233 13.40811892 12.69879233 13.40811892]
[12.70564858 13.41808752 12.70564858 13.41808752]
[12.71249166 13.42800632 12.71249166 13.42800632]
[12.71932138 13.43787576 12.71932138 13.43787576]
[12.72613754 13.44769627 12.72613754 13.44769627]
[12.73293995 13.45746827 12.73293995 13.45746827]
[12.73972843 13.46719218 12.73972843 13.46719218]
[12.74650282 13.47686839 12.74650282 13.47686839]
[12.75326294 13.4864973 12.75326294 13.4864973 ]
[12.76000864 13.49607931 12.76000864 13.49607931]
[12.76673977 13.50561479 12.76673977 13.50561479]
[12.77345618 13.51510411 12.77345618 13.51510411]
[12.78015773 13.52454765 12.78015773 13.52454765]
[12.78684429 13.53394577 12.78684429 13.53394577]
[12.79351573 13.54329881 12.79351573 13.54329881]
[12.80017193 13.55260713 12.80017193 13.55260713]
[12.80681277 13.56187107 12.80681277 13.56187107]
[12.81343814 13.57109096 12.81343814 13.57109096]
[12.82004793 13.58026713 12.82004793 13.58026713]]
plt.plot(tme_train, y_train, label = 'Training Data')
plt.plot(tme_test, np.exp(fcast_mean), label = 'Forecast', color = 'black') #Note the exponentiation
plt.plot(tme_test, y_test, color = 'red', label = 'Actual future values')
plt.plot(tme_test, np.exp(predlower), color = 'green', label = 'Prediction lower bound')
plt.plot(tme_test, np.exp(predupper), color = 'green', label = 'Prediction upper bound')
plt.axvline(x=n_train, color='gray', linestyle='--')
plt.legend()
plt.show()
Residuals: raw data versus logarithms¶
From the residual plots below, it is easy to see that the raw residuals grow more variable as house prices increase, while the log residuals have a much more stable spread. This supports roughly constant percentage errors rather than constant dollar errors. Compare how the spread changes over time, not the numerical sizes of residuals on these different scales.
resid_raw = AutoReg(y_train, lags=p).fit().resid
resid_log = armod_sm.resid
fig, axes = plt.subplots(2, 1, figsize=(11, 6), sharex=True)
for ax, resid, label in zip(axes, [resid_raw, resid_log],
['Raw data (dollars)', 'Log data (log units)']):
ax.plot(tme_train[p:], resid, color='black', linewidth=1)
ax.axhline(0, color='red', linestyle='--')
ax.set(title=label, ylabel='Residual')
axes[-1].set_xlabel('Quarter')
plt.tight_layout()
plt.show()
Residuals and AutoCorrelation¶
Given a set of values , how do we determine if is a good model for the data? Here are three example datasets.
n = 600
seed = 43
rng = np.random.default_rng(seed)
sig = 1
y_wn = rng.normal(loc = 0, scale = sig, size = n)
y_0 = np.concatenate(([0], y_wn))
fig, axes = plt.subplots(nrows = 3, ncols = 1, figsize = (12, 6))
axes[0].plot(y_wn)
axes[0].set_title('Dataset 1')
theta = 0.8
y_ma_1 = y_0[1:] + theta * y_0[:-1]
axes[1].plot(y_ma_1)
axes[1].set_title(f'Dataset 2')
theta = -0.8
y_ma_2 = y_0[1:] + theta * y_0[:-1] #y_ma_2 is simulated data from MA(1) with negative theta
axes[2].plot(y_ma_2)
axes[2].set_title(f'Dataset 3')
plt.tight_layout()
plt.show()
We can distinguish between these datasets by computing the sample ACF which is defined as follows. Given a time series dataset and a lag h, define
for . Here is simply the mean of . The acf plot graphs on the x-axis and on the y-axis. The quantity is known as the sample autocorrelation of the data at lag (acf stands for “Autocorrelation Function”).
Below we compute the sample ACF.
def sample_acf(y, h_max):
n = len(y)
y_mean = sum(y) / n
denominator = sum((y_t - y_mean) ** 2 for y_t in y)
autocorr = []
for h in range(h_max + 1):
numerator = sum((y[t] - y_mean) * (y[t + h] - y_mean) for t in range(n - h))
r_h = numerator / denominator
autocorr.append(r_h)
return autocorrThe sample ACF for each of the three datasets is plotted below.
fig, axes = plt.subplots(nrows = 3, ncols = 1, figsize = (12, 6))
h_max = 50
sample_acf_vals_wn = sample_acf(y_wn, h_max)
markerline, stemline, baseline = axes[0].stem(range(0, h_max + 1), sample_acf_vals_wn)
markerline.set_marker("None")
axes[0].set_title('Sample ACF of Dataset 1')
sample_acf_vals_ma_1 = sample_acf(y_ma_1, h_max)
markerline, stemline, baseline = axes[1].stem(range(0, h_max + 1), sample_acf_vals_ma_1)
markerline.set_marker("None")
theta = 0.8
axes[1].set_title(f'Simulated ACF of Dataset 2')
sample_acf_vals_ma_2 = sample_acf(y_ma_2, h_max)
markerline, stemline, baseline = axes[2].stem(range(0, h_max + 1), sample_acf_vals_ma_2)
markerline.set_marker("None")
theta = -0.8
axes[2].set_title(f'Simulated ACF of Dataset 3')
plt.tight_layout()
plt.show()
Note that the autocorrelation at lag 0 is always one. For the first dataset, all the other autocorrelations are negligible which suggests that this dataset is i.i.d N(0, ) (i.e., it is Gaussian white noise). For the other two datasets, the autocorrelation at lag 1 is nonnegligible. Therefore, the second and third datasets are not i.i.d .
There is an inbuilt function in statsmodels for plotting the sample_acf (and another inbuilt function for plotting the sample_acf). Its use is demonstrated below. Let us check that it gives identical answers to our computation.
from statsmodels.tsa.stattools import acf
acf_values_wn = acf(y_wn, nlags=h_max)
print(np.column_stack([sample_acf_vals_wn, acf_values_wn]))[[ 1.00000000e+00 1.00000000e+00]
[-6.18986853e-02 -6.18986853e-02]
[-3.16190776e-02 -3.16190776e-02]
[ 3.37684035e-02 3.37684035e-02]
[-5.08194511e-02 -5.08194511e-02]
[ 3.25331941e-02 3.25331941e-02]
[ 8.74858962e-03 8.74858962e-03]
[-1.16498128e-02 -1.16498128e-02]
[ 3.26472220e-02 3.26472220e-02]
[ 8.33178480e-02 8.33178480e-02]
[-3.27430872e-02 -3.27430872e-02]
[ 1.14743957e-02 1.14743957e-02]
[-6.46748355e-02 -6.46748355e-02]
[-9.26355222e-03 -9.26355222e-03]
[-1.48030902e-03 -1.48030902e-03]
[ 5.98188096e-03 5.98188096e-03]
[-4.20248527e-02 -4.20248527e-02]
[ 5.47736551e-02 5.47736551e-02]
[ 9.57612253e-03 9.57612253e-03]
[-5.55307936e-02 -5.55307936e-02]
[ 7.32655261e-02 7.32655261e-02]
[ 8.66336997e-03 8.66336997e-03]
[ 5.21590566e-02 5.21590566e-02]
[-1.51121264e-02 -1.51121264e-02]
[-2.06604710e-02 -2.06604710e-02]
[-2.69178655e-02 -2.69178655e-02]
[ 3.49261541e-02 3.49261541e-02]
[ 3.86381603e-02 3.86381603e-02]
[-2.07370593e-02 -2.07370593e-02]
[ 3.56664414e-02 3.56664414e-02]
[ 3.73042789e-02 3.73042789e-02]
[-1.58643239e-03 -1.58643239e-03]
[-2.62720921e-02 -2.62720921e-02]
[-7.90709599e-03 -7.90709599e-03]
[ 1.72123791e-02 1.72123791e-02]
[ 1.64248945e-02 1.64248945e-02]
[-4.29481074e-02 -4.29481074e-02]
[-4.08484611e-02 -4.08484611e-02]
[ 2.76715612e-02 2.76715612e-02]
[-1.97488528e-02 -1.97488528e-02]
[ 2.31475734e-03 2.31475734e-03]
[-4.12288416e-02 -4.12288416e-02]
[-1.05493982e-03 -1.05493982e-03]
[-1.60649886e-02 -1.60649886e-02]
[ 8.55000861e-06 8.55000861e-06]
[-2.36297467e-02 -2.36297467e-02]
[ 6.36622276e-02 6.36622276e-02]
[ 3.89111470e-02 3.89111470e-02]
[ 1.70688163e-02 1.70688163e-02]
[ 5.94731208e-02 5.94731208e-02]
[-1.23706214e-02 -1.23706214e-02]]
The following function plots the sample_acf values.
from statsmodels.graphics.tsaplots import plot_acf
fig, axes = plt.subplots(nrows = 3, ncols = 1, figsize = (12, 6))
plot_acf(y_wn, lags = h_max, ax = axes[0])
axes[0].set_title('Sample ACF of Dataset 1')
plot_acf(y_ma_1, lags = h_max, ax = axes[1])
theta = 0.8
axes[1].set_title(f'Simulated ACF of Dataset 2')
plot_acf(y_ma_2, lags = h_max, ax = axes[2])
theta = -0.8
axes[2].set_title(f'Simulated ACF of Dataset 3')
plt.tight_layout()
plt.show()
The plot_acf function gives a blue shading to assess whether each autocorrelation is non-negligible (significant) or not. Values appearing inside the blue area or just outside of it can be ignored.
Let us go back to the GDP dataset.
gdp = pd.read_csv("GDP_29Oct2026.csv")
print(gdp.head())
y = gdp['GDP'].to_numpy()
plt.plot(y, color = 'black')
plt.xlabel('Quarter')
plt.ylabel('Billions of Dollars')
plt.title('Gross Domestic Product')
plt.show() observation_date GDP
0 1947-01-01 243.164
1 1947-04-01 245.968
2 1947-07-01 249.585
3 1947-10-01 259.745
4 1948-01-01 265.742
n = len(y)
tme = range(1, n+1)
n_test = 50
#n_test = 12
n_train = n - n_test
y_train = y[:n_train]
tme_train = tme[:n_train]
y_test = y[n_train:]
tme_test = tme[n_train:]ylog_train = np.log(y_train)
ylog_test = np.log(y_test)We fit AR(1) to log-GDP.
from statsmodels.tsa.ar_model import AutoReg
armod_sm = AutoReg(ylog_train, lags = p).fit()
print(armod_sm.summary()) AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 269
Model: AutoReg(1) Log Likelihood 830.433
Method: Conditional MLE S.D. of innovations 0.011
Date: Fri, 30 Oct 2026 AIC -1654.866
Time: 22:18:16 BIC -1644.093
Sample: 1 HQIC -1650.539
269
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0283 0.004 7.232 0.000 0.021 0.036
y.L1 0.9984 0.000 2014.926 0.000 0.997 0.999
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.0016 +0.0000j 1.0016 0.0000
-----------------------------------------------------------------------------
To check whether this is a good model, we can plot the residuals.
resid_log = armod_sm.resid
fig, ax = plt.subplots(figsize=(11, 3))
ax.plot(tme_train[p:], resid_log, color='black', linewidth=1)
ax.axhline(0, color='red', linestyle='--')
ax.set(title='Log data (log units)', xlabel='Quarter', ylabel='Residual')
plt.tight_layout()
plt.show()
To assess whether the i.i.d assumption for the autoregression errors is justifiable, let us compute the sample acf of the residuals.
h_max = 50
fig, ax = plt.subplots()
plot_acf(resid_log, lags = h_max, ax = ax)
ax.set_title("Sample ACF of Residuals")
plt.show()
There are two non-negligible autocorrelations at lags 1 and 2. This suggests that i.i.d is not a good model for the residuals here. The Moving Average (MA) model of order 2 is more suitable for the residuals. We shall apply this model in the next lecture.