AutoRegressive (AR) Models

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

We will fit AutoRegressive Models to some time series datasets, and obtain predictions for future values of the time series.

Dataset One: Liquor Sales Data from FRED¶

#The following is FRED data on retail sales (in millions of dollars) for beer, wine and liquor stores (https://fred.stlouisfed.org/series/MRTSSM4453USN)
beersales = pd.read_csv('MRTSSM4453USN_October2025.csv')
print(beersales.head())
y = beersales['MRTSSM4453USN'].to_numpy()
plt.figure(figsize = (12, 6))
plt.plot(y)
plt.xlabel('Year')
plt.ylabel('Millions of Dollars')
plt.title('Retail Sales: Beer, wine and liquor stores')
plt.show()

  observation_date  MRTSSM4453USN
0       1992-01-01           1414
1       1992-02-01           1444
2       1992-03-01           1496
3       1992-04-01           1569
4       1992-05-01           1707

AutoRegressive (AR) models are a simple way for obtaining forecasts (predictions) for future values of the time series.

n = len(y)
print(n)
print(y)

403
[1414 1444 1496 1569 1707 1663 1792 1744 1658 1763 1715 2354 1512 1433
 1573 1606 1683 1679 1828 1665 1600 1647 1670 2287 1450 1410 1607 1648
 1715 1741 1858 1742 1710 1693 1731 2416 1459 1407 1619 1594 1696 1747
 1769 1745 1731 1678 1762 2437 1575 1549 1723 1687 1836 1837 1908 1936
 1672 1748 1871 2395 1605 1526 1745 1712 1940 1886 1976 1960 1783 1932
 1962 2657 1740 1665 1783 1848 2015 1938 2079 1993 1910 2033 2026 2857
 1755 1700 1867 1974 2078 2017 2200 2013 2000 2088 2140 3076 1810 1863
 2066 2011 2225 2278 2326 2288 2228 2223 2409 3192 1995 1941 2201 2110
 2337 2389 2349 2382 2191 2279 2495 3349 1999 1994 2236 2174 2429 2349
 2427 2451 2170 2267 2456 3294 2003 1931 2145 2202 2431 2297 2466 2511
 2323 2495 2514 3456 2169 2101 2270 2418 2543 2505 2716 2485 2481 2586
 2619 3604 2104 2194 2398 2501 2568 2637 2771 2645 2639 2656 2782 3925
 2303 2394 2616 2653 2870 2906 2966 2878 2852 2773 2987 4067 2445 2478
 2862 2730 3107 3203 3163 3136 2911 2942 3142 4195 2551 2677 2852 2861
 3264 3132 3359 3316 3010 3201 3235 4260 2782 2657 2888 2990 3312 3160
 3397 3249 3119 3274 3209 4422 2754 2790 3076 3168 3319 3292 3503 3307
 3222 3349 3374 4484 2765 2858 3111 3221 3325 3402 3584 3428 3373 3372
 3499 4696 2877 3077 3413 3278 3640 3668 3589 3606 3365 3482 3714 4917
 3021 3088 3465 3334 3771 3636 3848 3908 3459 3605 3863 4961 3213 3143
 3466 3511 3947 3822 4014 3988 3601 3857 3865 5245 3359 3273 3625 3631
 4070 3957 4252 3995 3818 4004 3995 5444 3376 3528 3786 3844 4114 4179
 4436 4182 4136 4138 4359 5645 3474 3503 3957 3960 4328 4353 4466 4324
 4221 4218 4506 5801 3684 3686 4314 4028 4574 4611 4632 4569 4296 4432
 4797 5969 3837 3818 4331 4252 4741 4647 4876 4902 4349 4569 4933 6163
 4066 4184 5081 4757 5737 5572 5905 5634 5410 5590 5631 7142 4960 4831
 5566 5560 5940 5877 6044 5601 5468 5647 5893 7474 4812 4920 5441 5544
 5727 5833 6003 5755 5621 5744 6017 7651 4983 4972 5535 5527 6067 6040
 6096 5834 5783 5736 6166 7863 4907 5248 5765 5539 6264 6155 6256 6203
 5658 5987 6453 7854 5055 5037 5613 5626 6310 5950 6225]

AR models work just like regression. The response vector and design matrix in the regression are set up as follows.

p = 24 #this is the order of the autoregression
yreg = y[p:] #these are the response values in the autoregression
Xmat = np.ones((n-p, 1)) #this will be the design matrix (X) in the autoregression
for j in range(1, p+1):
    col = y[p-j : n-j]
    Xmat = np.column_stack([Xmat, col])
print(Xmat.shape)
print(Xmat)

(379, 25)
[[1.000e+00 2.287e+03 1.670e+03 ... 1.496e+03 1.444e+03 1.414e+03]
 [1.000e+00 1.450e+03 2.287e+03 ... 1.569e+03 1.496e+03 1.444e+03]
 [1.000e+00 1.410e+03 1.450e+03 ... 1.707e+03 1.569e+03 1.496e+03]
 ...
 [1.000e+00 5.626e+03 5.613e+03 ... 6.096e+03 6.040e+03 6.067e+03]
 [1.000e+00 6.310e+03 5.626e+03 ... 5.834e+03 6.096e+03 6.040e+03]
 [1.000e+00 5.950e+03 6.310e+03 ... 5.783e+03 5.834e+03 6.096e+03]]

armod = sm.OLS(yreg, Xmat).fit()
print(armod.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.995
Model:                            OLS   Adj. R-squared:                  0.995
Method:                 Least Squares   F-statistic:                     3084.
Date:                Tue, 20 Oct 2026   Prob (F-statistic):               0.00
Time:                        19:00:07   Log-Likelihood:                -2274.2
No. Observations:                 379   AIC:                             4598.
Df Residuals:                     354   BIC:                             4697.
Df Model:                          24                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         10.1179     14.426      0.701      0.484     -18.254      38.490
x1             0.2560      0.052      4.932      0.000       0.154       0.358
x2             0.3601      0.054      6.706      0.000       0.255       0.466
x3             0.4011      0.057      7.038      0.000       0.289       0.513
x4            -0.0614      0.059     -1.032      0.303      -0.178       0.056
x5            -0.0263      0.060     -0.440      0.660      -0.144       0.091
x6            -0.0744      0.060     -1.230      0.219      -0.193       0.045
x7            -0.0822      0.060     -1.359      0.175      -0.201       0.037
x8            -0.0001      0.060     -0.002      0.998      -0.119       0.119
x9             0.2670      0.061      4.397      0.000       0.148       0.386
x10           -0.0916      0.059     -1.554      0.121      -0.208       0.024
x11            0.0157      0.056      0.279      0.781      -0.095       0.126
x12            0.7836      0.054     14.447      0.000       0.677       0.890
x13           -0.2662      0.054     -4.932      0.000      -0.372      -0.160
x14           -0.3699      0.056     -6.623      0.000      -0.480      -0.260
x15           -0.4112      0.059     -6.949      0.000      -0.528      -0.295
x16            0.0663      0.062      1.072      0.284      -0.055       0.188
x17            0.0471      0.062      0.760      0.448      -0.075       0.169
x18            0.0778      0.063      1.241      0.215      -0.046       0.201
x19            0.0813      0.063      1.296      0.196      -0.042       0.205
x20           -0.0005      0.063     -0.009      0.993      -0.124       0.123
x21           -0.2871      0.063     -4.550      0.000      -0.411      -0.163
x22            0.0834      0.061      1.360      0.175      -0.037       0.204
x23           -0.0047      0.058     -0.080      0.936      -0.119       0.110
x24            0.2444      0.056      4.348      0.000       0.134       0.355
==============================================================================
Omnibus:                       41.947   Durbin-Watson:                   1.964
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              131.720
Skew:                           0.463   Prob(JB):                     2.50e-29
Kurtosis:                       5.736   Cond. No.                     4.85e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.85e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

arfitted = armod.fittedvalues
#Plotting the fitted values against response values
plt.xlabel('Fitted Values')
plt.ylabel('Response Values')
plt.title('AR(1) Model Fitted Values vs Response Values')
plt.plot(arfitted, yreg, 'o')
plt.show()

#Generate k-step ahead forecasts: 
k = 100
yhat = np.concatenate([y, np.full(k, -9999)]) #extend data by k placeholder values
for i in range(1, k+1):
    ans = armod.params[0]
    for j in range(1, p+1):
        ans += armod.params[j] * yhat[n+i-j-1]
    yhat[n+i-1] = ans
predvalues = yhat[n:]

#Plotting the series with forecasts: 
plt.figure(figsize=(12, 6))
time_all = np.arange(1, n + k + 1)
plt.plot(time_all, yhat,  color='C0')
plt.plot(range(1, n + 1), y, label='Original Data', color='C1')
plt.plot(range(n + 1, n + k + 1), predvalues, label='Forecasts', color='blue')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

The predictions depend crucially on the order $p$ of the model. If $p \geq 12$ , we get reasonable predictions. However for smaller values of $p$ , the predictions look very unnatural. Go back and repeat the code above for other values of $p$ .

Dataset Two: House Price Data from FRED¶

#The following is FRED data on Average Sales Price of Houses Sold for the US (https://fred.stlouisfed.org/series/ASPUS)
hprice = pd.read_csv('ASPUS_October2025.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

p = 30 #this is the order of the AR model
n = len(y)
print(n)

yreg = y[p:] #these are the response values in the autoregression
Xmat = np.ones((n-p, 1)) #this will be the design matrix (X) in the autoregression
for j in range(1, p+1):
    col = y[p-j : n-j].reshape(-1, 1)
    Xmat = np.column_stack([Xmat, col])
print(Xmat.shape)
print(y)

(220, 31)
[ 19300  19400  19200  19600  19600  20200  20500  20900  21500  21000
  21600  21700  22700  23200  23000  22800  24000  24200  23900  24400
  25400  26700  26600  27000  27600  28100  28100  27100  27000  27300
  26000  26300  27300  28600  28300  28200  29000  29400  30300  31600
  32800  35100  35900  36600  38000  38600  39000  39300  40900  42600
  42200  44400  46000  47800  48100  50300  51600  54300  54000  57500
  59300  61600  63500  66400  68300  72400  74200  72700  73600  74400
  77500  80000  80900  84300  83800  83700  81200  85700  83900  84600
  86700  89100  92500  90800  94700  99200  98500  97800  98500 100500
 100500 103800 106300 112000 114400 115600 120800 126100 129900 133500
 137900 134800 141500 140400 144300 146800 150200 151200 149500 151200
 145500 150100 151100 148200 145400 144400 144500 145300 141700 147200
 144700 148900 148000 148300 153600 154200 152800 156100 153500 158900
 157700 160900 161100 166000 164000 171000 172200 177200 174700 175400
 180000 178800 184300 181500 189100 191800 193000 204800 202900 202400
 204100 212100 211000 211200 207800 214200 227600 227600 219100 232500
 233100 241000 248100 256000 262900 265300 274000 286300 288500 287800
 294600 294200 305300 302600 308100 299600 322100 310100 301200 305800
 290400 304200 285100 276600 257000 273400 274100 272900 275300 268800
 266000 278000 268100 267600 263000 259700 278000 282700 294500 297700
 307400 320400 324400 334400 331400 340600 340400 369400 348000 339700
 347400 366700 357000 357900 358800 364900 374800 376900 373200 399700
 374600 378400 392900 384000 375500 376700 382700 384600 383000 371100
 400600 396900 417400 428600 468000 496700 499300 525100 520300 521000
 505300 503000 521900 498300 519700 502200 498700 510900 514200 512800]

armod = sm.OLS(yreg, Xmat).fit()
print(armod.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.998
Model:                            OLS   Adj. R-squared:                  0.997
Method:                 Least Squares   F-statistic:                     2885.
Date:                Tue, 20 Oct 2026   Prob (F-statistic):          4.29e-235
Time:                        19:02:13   Log-Likelihood:                -2239.4
No. Observations:                 220   AIC:                             4541.
Df Residuals:                     189   BIC:                             4646.
Df Model:                          30                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       1252.3043    895.312      1.399      0.164    -513.784    3018.393
x1             0.7836      0.073     10.799      0.000       0.640       0.927
x2             0.5246      0.093      5.623      0.000       0.341       0.709
x3            -0.2076      0.100     -2.077      0.039      -0.405      -0.010
x4             0.0044      0.099      0.045      0.964      -0.191       0.200
x5             0.0292      0.099      0.296      0.768      -0.166       0.224
x6            -0.3315      0.098     -3.382      0.001      -0.525      -0.138
x7            -0.0284      0.101     -0.282      0.778      -0.227       0.170
x8             0.1678      0.102      1.646      0.101      -0.033       0.369
x9             0.1665      0.103      1.621      0.107      -0.036       0.369
x10           -0.2485      0.103     -2.413      0.017      -0.452      -0.045
x11            0.0576      0.106      0.546      0.586      -0.151       0.266
x12            0.3131      0.104      3.004      0.003       0.107       0.519
x13           -0.4220      0.108     -3.910      0.000      -0.635      -0.209
x14           -0.1454      0.110     -1.316      0.190      -0.363       0.073
x15            0.5566      0.116      4.807      0.000       0.328       0.785
x16           -0.0623      0.118     -0.527      0.599      -0.295       0.171
x17           -0.3085      0.123     -2.515      0.013      -0.550      -0.066
x18            0.1707      0.119      1.432      0.154      -0.064       0.406
x19            0.0860      0.118      0.728      0.468      -0.147       0.319
x20           -0.0764      0.118     -0.649      0.517      -0.309       0.156
x21           -0.2422      0.119     -2.043      0.042      -0.476      -0.008
x22            0.1030      0.119      0.867      0.387      -0.131       0.337
x23            0.0116      0.120      0.097      0.923      -0.225       0.248
x24           -0.2901      0.118     -2.468      0.014      -0.522      -0.058
x25            0.5416      0.117      4.641      0.000       0.311       0.772
x26            0.1503      0.123      1.219      0.224      -0.093       0.394
x27           -0.5970      0.123     -4.840      0.000      -0.840      -0.354
x28            0.2527      0.131      1.930      0.055      -0.006       0.511
x29            0.1942      0.121      1.610      0.109      -0.044       0.432
x30           -0.1401      0.103     -1.365      0.174      -0.343       0.062
==============================================================================
Omnibus:                       33.340   Durbin-Watson:                   1.992
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              114.033
Skew:                           0.550   Prob(JB):                     1.73e-25
Kurtosis:                       6.351   Cond. No.                     2.29e+06
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.29e+06. This might indicate that there are
strong multicollinearity or other numerical problems.

arfitted = armod.fittedvalues
#Plotting the fitted values against response values
plt.xlabel('Fitted Values')
plt.ylabel('Response Values')
plt.title('AR(1) Model Fitted Values vs Response Values')
#plotting with reduced size points for better visibility    
plt.plot(arfitted, yreg, 'o', markersize=4)  

plt.show()

#Generate k-step ahead forecasts: 
k = 100
yhat = np.concatenate([y, np.full(k, -9999)]) #extend data by k placeholder values
for i in range(1, k+1):
    ans = armod.params[0]
    for j in range(1, p+1):
        ans += armod.params[j] * yhat[n+i-j-1]
    yhat[n+i-1] = ans
predvalues = yhat[n:]

#Plotting the series with forecasts: 
time_all = np.arange(1, n + k + 1)
plt.plot(time_all, yhat, color='C0')
plt.plot(range(1, n + 1), y, label='Original Data', color='C1')
plt.plot(range(n + 1, n + k + 1), predvalues, label='Forecasts', color='blue')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

Repeat the above code for different values of $p$ to see how the predictions change.

We will now apply AR models to logarithms of $y$ .

ylog = np.log(y)
plt.plot(ylog)
plt.xlabel('Year')
plt.ylabel('Log Dollars')
plt.title('Log of Average Sales Price of Houses Sold in the US')
plt.show()

p = 30 #this is the order of the AR model
n = len(ylog)
print(n)

yreg = ylog[p:] #these are the response values in the autoregression
Xmat = np.ones((n-p, 1)) #this will be the design matrix (X) in the autoregression
for j in range(1, p+1):
    col = ylog[p-j : n-j].reshape(-1, 1)
    Xmat = np.column_stack([Xmat, col])

armod_log = sm.OLS(yreg, Xmat).fit()
print(armod_log.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                     7003.
Date:                Tue, 20 Oct 2026   Prob (F-statistic):          1.96e-271
Time:                        19:03:28   Log-Likelihood:                 509.83
No. Observations:                 220   AIC:                            -957.7
Df Residuals:                     189   BIC:                            -852.5
Df Model:                          30                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.1013      0.033      3.049      0.003       0.036       0.167
x1             0.8024      0.073     11.057      0.000       0.659       0.946
x2             0.4279      0.092      4.628      0.000       0.246       0.610
x3            -0.1316      0.097     -1.356      0.177      -0.323       0.060
x4             0.0597      0.095      0.627      0.531      -0.128       0.247
x5            -0.0348      0.095     -0.366      0.715      -0.222       0.153
x6            -0.2380      0.094     -2.544      0.012      -0.423      -0.053
x7             0.0186      0.095      0.195      0.845      -0.169       0.207
x8             0.0834      0.095      0.874      0.383      -0.105       0.272
x9             0.0662      0.095      0.694      0.488      -0.122       0.254
x10           -0.1635      0.095     -1.727      0.086      -0.350       0.023
x11            0.0860      0.095      0.905      0.367      -0.101       0.273
x12            0.2016      0.094      2.135      0.034       0.015       0.388
x13           -0.2995      0.096     -3.128      0.002      -0.488      -0.111
x14           -0.0542      0.097     -0.558      0.577      -0.246       0.137
x15            0.3066      0.097      3.161      0.002       0.115       0.498
x16           -0.1182      0.097     -1.214      0.226      -0.310       0.074
x17           -0.1742      0.098     -1.781      0.077      -0.367       0.019
x18            0.0710      0.097      0.735      0.463      -0.119       0.261
x19            0.1175      0.096      1.227      0.221      -0.071       0.306
x20            0.1415      0.096      1.473      0.143      -0.048       0.331
x21           -0.2030      0.096     -2.105      0.037      -0.393      -0.013
x22            0.0006      0.097      0.007      0.995      -0.191       0.192
x23           -0.0754      0.097     -0.780      0.437      -0.266       0.115
x24           -0.1470      0.097     -1.517      0.131      -0.338       0.044
x25            0.3098      0.096      3.228      0.001       0.120       0.499
x26            0.1047      0.098      1.065      0.288      -0.089       0.299
x27           -0.3197      0.098     -3.256      0.001      -0.513      -0.126
x28            0.1514      0.101      1.503      0.135      -0.047       0.350
x29            0.0733      0.097      0.756      0.451      -0.118       0.264
x30           -0.0703      0.075     -0.938      0.349      -0.218       0.078
==============================================================================
Omnibus:                        0.073   Durbin-Watson:                   1.980
Prob(Omnibus):                  0.964   Jarque-Bera (JB):                0.002
Skew:                           0.001   Prob(JB):                        0.999
Kurtosis:                       3.013   Cond. No.                     6.47e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.47e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

arfitted_log = armod_log.fittedvalues
#Plotting the fitted values against response values
plt.xlabel('Fitted Values')
plt.ylabel('Response Values')
plt.title('AR(1) Model Fitted Values vs Response Values')
#plotting with reduced size points for better visibility    
plt.plot(arfitted_log, yreg, 'o', markersize=4)  

plt.show()

#Generate k-step ahead forecasts: 
k = 100
yhat = np.concatenate([ylog, np.full(k, -9999)]) #extend data by k placeholder values
for i in range(1, k+1):
    ans = armod_log.params[0]
    for j in range(1, p+1):
        ans += armod_log.params[j] * yhat[n+i-j-1]
    yhat[n+i-1] = ans
predvalues_log = yhat[n:]

#Plotting the series with forecasts: 
time_all = np.arange(1, n + k + 1)
plt.plot(time_all, yhat, color='C0')
plt.plot(range(1, n + 1), ylog, label='Original Data (log scale)', color='C1')
plt.plot(range(n + 1, n + k + 1), predvalues_log, label='Forecasts', color='blue')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

We now plot the predictions on the original data.

#Plotting the series with forecasts: 
time_all = np.arange(1, n + k + 1)
plt.plot(range(1, n + 1), y, label='Original Data', color='C1')
plt.plot(range(n + 1, n + k + 1), np.exp(predvalues_log), label='Forecasts (AR on log)', color='blue')
plt.plot(range(n + 1, n + k + 1), predvalues, label='Forecasts (AR original data)', color='green')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

Sunspots Data¶

Autoregressive models were invented in the context of the sunspots data by Yule in 1927.

sunspots = pd.read_csv('SN_y_tot_V2.0_17Sept2026.csv', header = None, sep = ';')
y = sunspots.iloc[:,1].values
n = len(y)
plt.figure(figsize = (12, 6))
plt.plot(y)
plt.show()

The most basic model for this dataset is the single sinusoid model:

y_t = \beta_0 + \beta_1 \cos(2 \pi f t) + \beta_2 \sin(2 \pi f t) + \epsilon_t

with $\epsilon_t \overset{\text{i.i.d}}{\sim} N(0, \sigma^2)$ . We previously used this model. For estimating $f$ , we calculate $RSS(f)$ for a bunch of $f$ values, and then minimize it over $f$ .

#rss function:
def rss(f):
    n = len(y)
    X = np.column_stack([np.ones(n)])
    x = np.arange(1, n+1)
    xcos = np.cos(2 * np.pi * f * x)
    xsin = np.sin(2 * np.pi * f * x)
    X = np.column_stack([X, xcos, xsin])
    md = sm.OLS(y, X).fit()
    ans = np.sum(md.resid ** 2)
    return ans

ngrid = 10000
fvals = np.linspace(0, 0.5, ngrid)
rssvals = np.array([rss(f) for f in fvals])
plt.plot(fvals, rssvals) #plot rss and find f which minimizes rss
plt.show()

fhat = fvals[np.argmin(rssvals)]
print(fhat) 
print(1/fhat)  #period corresponding to fhat

0.09090909090909091
11.0

The MLE of $f$ is basically $1/11$ (corresponding to the 11-year solar cycle). After estimating $f$ , the other parameters $\beta_0, \beta_1, \beta_2, \sigma$ are estimated as in usual linear regression as follows.

#Estimates of beta and sigma: 
x = np.arange(1, n+1)
f = fhat
xcos = np.cos(2 * np.pi * f * x)
xsin = np.sin(2 * np.pi * f * x)
X = np.column_stack([np.ones(n), xcos, xsin])
md = sm.OLS(y, X).fit()
print(md.params) #this gives estimates of beta_0, beta_1, beta_2 
rss_fhat = np.sum(md.resid ** 2)
#there are two estimates for sigma (which usually give similar values)
sigma_mle = np.sqrt(rss_fhat/n)
sigma_unbiased = np.sqrt((rss_fhat)/(n-3))
print(np.array([sigma_mle, sigma_unbiased]))

[ 78.86715227 -38.27095266 -29.26525055]
[51.54476939 51.78358813]

Back in Lecture 7, we also computed an uncertainty interval (credible region) for $f$ , using the following posterior formula:

\begin{align*} \text{posterior}(f) \propto I\{0 < f < 1/2 \} |X_f^T X_f|^{-1/2} \left(\frac{1}{RSS(f)} \right)^{(n-p)/2} \end{align*}

(1)

with $p = 3$ (the number of columns of $X_f$ ). We evaluated this posterior on a fine grid of frequencies as in the code below.

def logpost(f):
    x = np.arange(1, n+1)
    xcos = np.cos(2 * np.pi * f * x)
    xsin = np.sin(2 * np.pi * f * x)
    X = np.column_stack([np.ones(n), xcos, xsin])
    p = X.shape[1]
    md = sm.OLS(y, X).fit()
    rss = np.sum(md.resid ** 2)
    sgn, log_det = np.linalg.slogdet(np.dot(X.T, X)) #sgn gives the sign of the determinant (in our case, this should 1)
    #log_det gives the logarithm of the absolute value of the determinant
    logval = ((p-n)/2) * np.log(rss) - (0.5)*log_det
    return logval

allfvals = np.arange(0.01, 0.5, 0.00001) #fine grid of frequencies (staying away from 0 where the determinant term is close to zero)
logpostvals = np.array([logpost(f) for f in allfvals])
postvals = np.exp(logpostvals - np.max(logpostvals))
postvals = postvals/(np.sum(postvals))
print(allfvals[np.argmax(postvals)]) #posterior mode (same as the RSS minimizer)
print(1/allfvals[np.argmax(postvals)])

0.09088999999999671
11.00231048520229

#95% credible interval: smallest symmetric region (in the grid) around the posterior mode with posterior probability at least 0.95
def PostProbAroundMax(m):
    est_ind = np.argmax(postvals)
    ans = np.sum(postvals[(est_ind-m):(est_ind+m+1)])
    return(ans)
m = 0
while PostProbAroundMax(m) <= 0.95:
    m = m+1
est_ind = np.argmax(postvals)
f_est = allfvals[est_ind]
#95% credible interval for f:
ci_f_low = allfvals[est_ind - m]
ci_f_high = allfvals[est_ind + m]
print(np.array([f_est, ci_f_low, ci_f_high]))

#For the period (1/frequency):
period_est = 1/f_est
ci_period_low = 1/ci_f_high
ci_period_high = 1/ci_f_low
print(np.array([period_est, ci_period_low, ci_period_high]))
#uncertainty in days:
print(np.array([(period_est - ci_period_low)*365, (ci_period_high - period_est)*365]))

[0.09089 0.09061 0.09117]
[11.00231049 10.96852035 11.03630946]
[12.33340059 12.40962511]

As in Lecture 7, the 95% credible interval for the period is $[10.97, 11.04]$ years, i.e.,

\begin{align*} [11~\text{years} - 12~\text{days}, 11~\text{years} + 12~\text{days}]. \end{align*}

(2)

Note that this interval is quite narrow.

Yule did not like the simple single sinusoidal model for the sunspots data, mainly because if we simulate data from it, the data will look much more noisy and irregular as opposed to the sunspots data which have a smooth appearance. We can check this by simulating synthetic data from the single sinusoid model and comparing the synthetic data with the actual sunspots data. We will take the parameters $f, \beta, \sigma$ to be those estimated from the sunspots dataset to make the plots comparable to the actual sunspots data.

rng = np.random.default_rng(seed = 42)
errorsamples = rng.normal(loc = 0, scale = sigma_mle, size = n)
sim_data = md.fittedvalues + errorsamples
plt.figure(figsize = (10, 6))
plt.plot(sim_data)
plt.show()

To facilitate comparison with the actual sunspots dataset, let us plot a bunch of these simulated datasets along with the real sunspots data (just to see if the sunspots dataset can be spotted as the “odd one out” from these plots).

fig, axes = plt.subplots(3, 3, figsize = (12, 4))
axes = axes.flatten()
for i in range(5):
    errorsamples = rng.normal(loc = 0, scale = sigma_mle, size = n)
    sim_data = md.fittedvalues + errorsamples
    axes[i].plot(sim_data)
axes[5].plot(y)
for i, idx in enumerate(range(6, 9)):
    errorsamples = rng.normal(loc = 0, scale = sigma_mle, size = n)
    sim_data = md.fittedvalues + errorsamples
    axes[idx].plot(sim_data)
plt.tight_layout()
plt.show()

It is clear from the above plot-grid that the sunspots dataset sticks out as the odd one out. This shows that the simulated datasets (which are too wiggly and do not have well-defined peaks) are generated from a model that ignores many aspects of the sunspots dataset.

Yule Model¶

Yule model is given by:

y_t = \phi_0 + \phi_1 y_{t-1} - y_{t-2} + \epsilon_t

Note that $\phi_2$ in the usual AR(2) model is set to -1 above. Also $\phi_1$ is given by $2 \cos(2 \pi f)$ . The above model can be rewritten as

y_t + y_{t-2} = \phi_0 + \phi_1 y_{t-1} + \epsilon_t

so we can estimate $\phi_0$ and $\phi_1$ by regressing $y_t + y_{t-2}$ on $y_{t-1}$ .

p = 2
yreg = y[p:]
x1 = y[1:-1]
x2 = y[:-2]
Xmat = np.column_stack([np.ones(len(yreg)), x1])
print(Xmat.shape)

(324, 2)

y_adjusted = yreg + x2
yulemod = sm.OLS(y_adjusted, Xmat).fit()
print(yulemod.summary())
print(yulemod.params)
slpe_est = yulemod.params[1]

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.930
Model:                            OLS   Adj. R-squared:                  0.929
Method:                 Least Squares   F-statistic:                     4258.
Date:                Tue, 20 Oct 2026   Prob (F-statistic):          1.06e-187
Time:                        19:39:27   Log-Likelihood:                -1537.0
No. Observations:                 324   AIC:                             3078.
Df Residuals:                     322   BIC:                             3086.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         28.7410      2.513     11.439      0.000      23.798      33.684
x1             1.6345      0.025     65.254      0.000       1.585       1.684
==============================================================================
Omnibus:                       11.183   Durbin-Watson:                   2.484
Prob(Omnibus):                  0.004   Jarque-Bera (JB):               16.872
Skew:                           0.237   Prob(JB):                     0.000217
Kurtosis:                       4.013   Cond. No.                         163.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[28.74104583  1.63446369]

From the estimate for $\phi_1$ , the frequency parameter can be estimated using the relation: $\phi_1 = 2 \cos (2 \pi f)$ .

yulefhat = np.arccos(slpe_est/2) / (2 * np.pi)
print(1/yulefhat)

10.229784670602625

The period corresponding to this frequency is 10.259 which is smaller than the period of 11 that we get by fitting $y_t = \beta_0 + \beta_1 \cos(2 \pi f t) + \beta_2 \sin(2 \pi f t) + \epsilon_t$ .

Uncertainty quantification for $f$ in the Yule model¶

How uncertain is this estimate of the period? The Yule model is a linear regression (of $z_t = y_t + y_{t-2}$ on $y_{t-1}$ with an intercept), so we can use the Bayesian analysis of linear regression from Lecture 4. With the flat prior on $(\phi_0, \phi_1)$ and the prior $1/\sigma$ on $\sigma$ , the posterior is

\begin{align*} \sigma^2 \mid \text{data} \sim \frac{RSS}{\chi^2_{N - p}} \qquad \text{and} \qquad (\phi_0, \phi_1) \mid \sigma, \text{data} \sim N\left(\hat{\phi}, \sigma^2 (X^T X)^{-1}\right) \end{align*}

(3)

where $N = n - 2$ is the number of observations in the regression and $p = 2$ is the number of parameters. Integrating out $\sigma$ , the marginal posterior of $\phi_1$ is a shifted and scaled $t$ -distribution:

\begin{align*} \phi_1 \mid \text{data} \sim \hat{\phi}_1 + \hat{\sigma} \sqrt{(X^T X)^{-1}_{11}} \; t_{N-p} \qquad \text{with } \hat{\sigma}^2 = \frac{RSS}{N - p}. \end{align*}

(4)

This gives a 95% credible interval for $\phi_1$ . Since $f = \arccos(\phi_1/2)/(2 \pi)$ is a decreasing function of $\phi_1$ (and the period $1/f$ is an increasing function of $\phi_1$ ), the 95% credible interval for $f$ (or for the period) is obtained by simply transforming the endpoints of the interval for $\phi_1$ .

from scipy import stats

N_yule, p_yule = Xmat.shape #N = n - 2 observations and p = 2 parameters in the Yule regression
phi1_hat = yulemod.params[1]
rss_yule = np.sum(yulemod.resid ** 2)
sigma2_hat = rss_yule/(N_yule - p_yule)
XtX_inv = np.linalg.inv(Xmat.T @ Xmat)
se_phi1 = np.sqrt(sigma2_hat * XtX_inv[1, 1]) #this is the "std err" of x1 in the regression summary above
print(np.array([phi1_hat, se_phi1]))

#95% credible interval for phi1 (shifted and scaled t distribution with N - p degrees of freedom):
tq = stats.t.ppf(0.975, N_yule - p_yule)
phi1_low = phi1_hat - tq * se_phi1
phi1_high = phi1_hat + tq * se_phi1
print(np.array([phi1_hat, phi1_low, phi1_high])) #these coincide with the 95% confidence interval in the regression summary

[1.63446369 0.02504764]
[1.63446369 1.58518601 1.68374138]

#f = arccos(phi1/2)/(2 pi) is decreasing in phi1, so the endpoints get swapped:
yule_f_est = np.arccos(phi1_hat/2) / (2 * np.pi)
yule_ci_f_low = np.arccos(phi1_high/2) / (2 * np.pi)
yule_ci_f_high = np.arccos(phi1_low/2) / (2 * np.pi)
print(np.array([yule_f_est, yule_ci_f_low, yule_ci_f_high])) #95% credible interval for f in the Yule model

#For the period (1/frequency), as before:
yule_period_est = 1/yule_f_est
yule_ci_period_low = 1/yule_ci_f_high
yule_ci_period_high = 1/yule_ci_f_low
print(np.array([yule_period_est, yule_ci_period_low, yule_ci_period_high]))
#uncertainty in days:
print(np.array([(yule_period_est - yule_ci_period_low)*365, (yule_ci_period_high - yule_period_est)*365]))

[0.09775377 0.09072725 0.10436529]
[10.22978467  9.58172996 11.02204688]
[236.53996816 289.17570528]

So the 95% credible interval for the period under the Yule model is $[9.58, 11.02]$ years, i.e.,

\begin{align*} [10.23~\text{years} - 237~\text{days}, 10.23~\text{years} + 289~\text{days}]. \end{align*}

(5)

Compare this with the interval $[10.97, 11.04]$ years (11 years $\pm 12$ days) from the single sinusoid model: the Yule interval is about 40 times wider, and it (just) contains 11.

Why are the two intervals so different? In the single sinusoid model, $f$ is determined by phase coherence over the whole dataset: if $f$ is slightly off, the fitted sinusoid drifts out of phase with the data over the roughly 30 cycles in the dataset, and $RSS(f)$ increases sharply. In the Yule model, $\hat{\phi}_1$ is essentially read off the local relation between $y_t + y_{t-2}$ and $y_{t-1}$ , which is much noisier.

#Generate k-step ahead forecasts: 
k = 100
yhat = np.concatenate([y, np.full(k, -9999)]) #extend data by k placeholder values
for i in range(1, k+1):
    ans = yulemod.params[0] - yhat[n+i-3]
    ans += yulemod.params[1] * yhat[n+i-2]
    yhat[n+i-1] = ans
predvalues = yhat[n:]
print(predvalues)
print(yhat)

[ 75.40697305  28.79100562   0.39192621   0.59062938  29.31448189
  76.06387284 123.7502026  154.94238637 158.23854857 132.43382223
  86.96077169  38.4414478    4.61142496  -2.16319529  20.5939567
  64.56431569 113.67511911 149.97458535 160.19394164 140.59764224
  98.34884602  48.89102186  10.30280005  -3.3104234   13.02747893
  53.34441057 102.90306931 143.58796615 160.52729422 147.52911413
 109.3447326   59.93192737  17.35287269  -2.82824113   6.7655157
  42.62727674  91.64826639 135.9095332  159.23197725 153.09039853
 119.72976701  71.34460469  25.62144501  -0.72623717   1.93259253
  32.62603532  80.13452355 127.09197996 156.33374945 157.17090363
 129.29743227  82.9021011   34.94408803   2.95378798  -1.37518299
  23.53957119  68.59080332 117.31065247 151.89024501 159.68948444
 137.85746561  94.374584    45.1354115    8.13865328  -3.09203236
  15.54857792  57.24666431 106.76006238 145.98982755 160.59605641
 145.23964204 105.53391139  55.99275055  14.72525239  -3.18381428
   8.81196458  46.3276963   95.65001893 138.74983289 159.87259143
 151.29715946 116.15816868  67.30019595  22.58260409  -1.64870359
   3.46369557  36.05103409  84.20145664 130.31423568 157.53347635
 155.90955797 126.03608169  78.83288763  31.55445694   1.48277248
  -0.38987332  26.62103957  72.64204184 120.85078638 153.62522683]
[  8.3         18.3         26.7         38.3         60.
  96.7         48.3         33.3         16.7         13.3
   5.           0.           0.           3.3         18.3
  45.          78.3        105.         100.          65.
  46.7         43.3         36.7         18.3         35.
  66.7        130.         203.3        171.7        121.7
  78.3         58.3         18.3          8.3         26.7
  56.7        116.7        135.         185.         168.3
 121.7         66.7         33.3         26.7          8.3
  18.3         36.7         66.7        100.         134.8
 139.          79.5         79.7         51.2         20.3
  16.          17.          54.          79.3         90.
 104.8        143.2        102.          75.2         60.7
  34.8         19.          63.         116.3        176.8
 168.         136.         110.8         58.          51.
  11.7         33.         154.2        257.3        209.8
 141.3        113.5         64.2         38.          17.
  40.2        138.2        220.         218.2        196.8
 149.8        111.         100.          78.2         68.3
  35.5         26.7         10.7          6.8         11.3
  24.2         56.7         75.          71.8         79.2
  70.3         46.8         16.8         13.5          4.2
   0.           2.3          8.3         20.3         23.2
  59.          76.3         68.3         52.9         38.5
  24.2          9.2          6.3          2.2         11.4
  28.2         59.9         83.         108.5        115.2
 117.4         80.8         44.3         13.4         19.5
  85.8        192.7        227.3        168.7        143.
 105.5         63.3         40.3         18.1         25.1
  65.8        102.7        166.3        208.3        182.5
 126.3        122.         102.7         74.1         39.
  12.7          8.2         43.4        104.4        178.3
 182.2        146.6        112.1         83.5         89.2
  57.8         30.7         13.9         62.8        123.6
 232.         185.3        169.2        110.1         74.5
  28.3         18.9         20.7          5.7         10.
  53.7         90.5         99.         106.1        105.8
  86.3         42.4         21.8         11.2         10.4
  11.8         59.5        121.7        142.         130.
 106.6         69.4         43.8         44.4         20.2
  15.7          4.6          8.5         40.8         70.1
 105.5         90.1        102.8         80.9         73.2
  30.9          9.5          6.           2.4         16.1
  79.          95.         173.6        134.6        105.7
  62.7         43.5         23.7          9.7         27.9
  74.         106.5        114.7        129.7        108.2
  59.4         35.1         18.6          9.2         14.6
  60.2        132.8        190.6        182.6        148.
 113.          79.2         50.8         27.1         16.1
  55.3        154.3        214.7        193.         190.7
 118.9         98.3         45.          20.1          6.6
  54.2        200.7        269.3        261.7        225.1
 159.          76.4         53.4         39.9         15.
  22.          66.8        132.9        150.         149.4
 148.          94.4         97.6         54.1         49.2
  22.5         18.4         39.3        131.         220.1
 218.9        198.9        162.4         91.          60.5
  20.6         14.8         33.9        123.         211.1
 191.8        203.3        133.          76.1         44.9
  25.1         11.6         28.9         88.3        136.3
 173.9        170.4        163.6         99.3         65.3
  45.8         24.7         12.6          4.2          4.8
  24.9         80.8         84.5         94.         113.3
  69.8         39.8         21.7          7.           3.6
   8.8         29.6         83.2        125.5        154.7
 123.2         75.40697305  28.79100562   0.39192621   0.59062938
  29.31448189  76.06387284 123.7502026  154.94238637 158.23854857
 132.43382223  86.96077169  38.4414478    4.61142496  -2.16319529
  20.5939567   64.56431569 113.67511911 149.97458535 160.19394164
 140.59764224  98.34884602  48.89102186  10.30280005  -3.3104234
  13.02747893  53.34441057 102.90306931 143.58796615 160.52729422
 147.52911413 109.3447326   59.93192737  17.35287269  -2.82824113
   6.7655157   42.62727674  91.64826639 135.9095332  159.23197725
 153.09039853 119.72976701  71.34460469  25.62144501  -0.72623717
   1.93259253  32.62603532  80.13452355 127.09197996 156.33374945
 157.17090363 129.29743227  82.9021011   34.94408803   2.95378798
  -1.37518299  23.53957119  68.59080332 117.31065247 151.89024501
 159.68948444 137.85746561  94.374584    45.1354115    8.13865328
  -3.09203236  15.54857792  57.24666431 106.76006238 145.98982755
 160.59605641 145.23964204 105.53391139  55.99275055  14.72525239
  -3.18381428   8.81196458  46.3276963   95.65001893 138.74983289
 159.87259143 151.29715946 116.15816868  67.30019595  22.58260409
  -1.64870359   3.46369557  36.05103409  84.20145664 130.31423568
 157.53347635 155.90955797 126.03608169  78.83288763  31.55445694
   1.48277248  -0.38987332  26.62103957  72.64204184 120.85078638
 153.62522683]

#Plotting the series with forecasts: 
plt.figure(figsize=(12, 6))
time_all = np.arange(1, n + k + 1)
plt.plot(time_all, yhat, label='Extended Series (yhat)', color='C0')
plt.plot(range(1, n + 1), y, label='Original Data', color='C1')
plt.plot(range(n + 1, n + k + 1), predvalues, label='Forecasts', color='blue')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

The predictions above are perfectly sinusoidal. Because the cycles for the sunspots data are irregular and can have different periods, there is a danger that these perfectly sinusoidal predictions can get out of phase sometimes leading to loss of prediction accuracy.

Below, we simulate datasets using the Yule model.

#Simulating from the Yule Model:
rng = np.random.default_rng(seed = 42)
sighat = np.sqrt(np.mean(yulemod.resid ** 2))
print(sighat)
fig, axes = plt.subplots(3, 3, figsize = (12, 6))
axes = axes.flatten()
for j in range(1):
    ysim = y.copy()
    for i in range(2, n):
        err = rng.normal(loc = 0, scale = sighat, size = 1)
        ysim[i] = yulemod.params[0] + yulemod.params[1] * ysim[i-1] - ysim[i-2] + err[0]
    axes[j].plot(ysim)
axes[1].plot(y)    
for j in range(2, 9):
    ysim = y.copy()
    for i in range(2, n):
        err = rng.normal(loc = 0, scale = sighat, size = 1)
        ysim[i] = yulemod.params[0] + yulemod.params[1] * ysim[i-1] - ysim[i-2] + err[0]
    axes[j].plot(ysim)
plt.show()

27.796017406459733

The sunspots dataset still probably sticks out compared to the synthetic datasets. However, all these synthetic datasets are quite smooth just like the actual sunspots data. This plots shows that the two “sinusoid + noise” models described in lecture are fundamentally different.

AR(2) Model¶

The Yule model is a special case of the AR(2) model. Specifically, the AR(2) model is: $y_t = \phi_0 + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \epsilon_t$ and the Yule model is a special case obtained by taking $\phi_2 = 1$ . Below we fit the AR(2) model to the sunspots data (i.e., we estimate all the parameters $\phi_0, \phi_1, \phi_2$ ).

p = 2 #this is the order of the AR model
n = len(y)
print(n)

yreg = y[p:] #these are the response values in the autoregression
Xmat = np.ones((n-p, 1)) #this will be the design matrix (X) in the autoregression
for j in range(1, p+1):
    col = y[p-j : n-j].reshape(-1, 1)
    Xmat = np.column_stack([Xmat, col])

armod = sm.OLS(yreg, Xmat).fit()
print(armod.params)
print(armod.summary())
sighat = np.sqrt(np.mean(armod.resid ** 2))
print(sighat)

[24.49130086  1.38543496 -0.69541909]
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.828
Model:                            OLS   Adj. R-squared:                  0.827
Method:                 Least Squares   F-statistic:                     775.0
Date:                Tue, 20 Oct 2026   Prob (F-statistic):          1.34e-123
Time:                        19:44:10   Log-Likelihood:                -1510.3
No. Observations:                 324   AIC:                             3027.
Df Residuals:                     321   BIC:                             3038.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         24.4913      2.384     10.272      0.000      19.801      29.182
x1             1.3854      0.040     34.506      0.000       1.306       1.464
x2            -0.6954      0.040    -17.315      0.000      -0.774      -0.616
==============================================================================
Omnibus:                       33.598   Durbin-Watson:                   2.195
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               50.854
Skew:                           0.671   Prob(JB):                     9.06e-12
Kurtosis:                       4.401   Cond. No.                         232.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
25.59744475509888

#Simulations from the fitted AR(2) model:
fig, axes = plt.subplots(3, 3, figsize = (12, 6))
axes = axes.flatten()

for j in range(3):
    ysim = y.copy()
    for i in range(2, n):
        err = rng.normal(loc = 0, scale = sighat, size = 1)
        ysim[i] = armod.params[0] + armod.params[1] * ysim[i-1] + armod.params[2] * ysim[i-2] + err[0]
    axes[j].plot(ysim)

axes[3].plot(y)

for j in range(4, 9):
    ysim = y.copy()
    for i in range(2, n):
        err = rng.normal(loc = 0, scale = sighat, size = 1)
        ysim[i] = armod.params[0] + armod.params[1] * ysim[i-1] + armod.params[2] * ysim[i-2] + err[0]
    axes[j].plot(ysim)
plt.show()

These datasets are still smooth (much smoother than observations simulated from $y_t = \beta_0 + \beta_1 \cos(2 \pi f t) + \beta_2 \sin(2 \pi f t) + \epsilon_t$ ).

Below, we obtain predictions using the AR(2) model.

#Generate k-step ahead forecasts: 
k = 100
yhat = np.concatenate([y, np.full(k, -9999)]) #extend data by k placeholder values
for i in range(1, k+1):
    ans = armod.params[0]
    for j in range(1, p+1):
        ans += armod.params[j] * yhat[n+i-j-1]
    yhat[n+i-1] = ans
predvalues = yhat[n:]

#Plotting the series with forecasts: 
plt.figure(figsize=(12, 6))
time_all = np.arange(1, n + k + 1)
plt.plot(time_all, yhat, color='C0')
plt.plot(range(1, n + 1), y, label='Original Data', color='C1')
plt.plot(range(n + 1, n + k + 1), predvalues, label='Forecasts', color='blue')
plt.axvline(x=n, color='black', linestyle='--', label='Forecast Start')
#plt.axhline(y=np.mean(y), color='gray', linestyle=':', label='Mean of Original Data')
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('Time Series + AR(p) Forecasts')
plt.legend()
plt.show()

The predictions above are quite different from those obtained from the Yule model. The predictions actually are given by a damped sinusoid (as we shall see in the coming lectures). The predictions become flat after a few time points. When the cycles become irregular with varying periods, constant predictors might do well compared to sinusoidal predictors which may go out of phase at some places.

Here are all the three predictions in one plot.

k = 100
future_t = np.arange(n + 1, n + k + 1)

# Sinusoid plus noise: predict the fitted mean.
pred_sinusoid = (md.params[0]
                 + md.params[1] * np.cos(2 * np.pi * fhat * future_t)
                 + md.params[2] * np.sin(2 * np.pi * fhat * future_t))

# Yule and AR(2): predict recursively, setting future errors to zero.
pred_yule = np.concatenate([y, np.zeros(k)])
pred_ar2 = np.concatenate([y, np.zeros(k)])
for t in range(n, n + k):
    pred_yule[t] = yulemod.params[0] + yulemod.params[1] * pred_yule[t-1] - pred_yule[t-2]
    pred_ar2[t] = armod.params[0] + armod.params[1] * pred_ar2[t-1] + armod.params[2] * pred_ar2[t-2]

plt.figure(figsize=(14, 6))
plt.plot(np.arange(1, n + 1), y, color='black', label='Sunspots data')
plt.plot(future_t, pred_sinusoid, color='blue', label='Sinusoid plus noise')
plt.plot(future_t, pred_yule[n:], color='red', label='Yule model')
plt.plot(future_t, pred_ar2[n:], color='green', label='AR(2)')
plt.axvline(n, color='gray', linestyle='--')
plt.xlabel('Time (years)')
plt.ylabel('Sunspot number')
plt.title('Sunspots data and forecasts for the next 100 years')
plt.legend()
plt.tight_layout()
plt.show()

CPI Data¶

We can obtain the same distinction (of placing the noise in two different places) in other contexts as well. Below, we illustrate this phenomenon in the context of a linear trend model. We will look at the example of CPI data.

We shall use the Consumer Price Index (CPI) Data (downloaded from https://fred.stlouisfed.org/series/CPIAUCSL). Changes in the CPI are used to measure Inflation.

## Inflation (Consumer Price Index)
cpi = pd.read_csv('CPIAUCSL_01September2025.csv')
print(cpi)

    observation_date  CPIAUCSL
0         1947-01-01    21.480
1         1947-02-01    21.620
2         1947-03-01    22.000
3         1947-04-01    22.000
4         1947-05-01    21.950
..               ...       ...
938       2025-03-01   319.615
939       2025-04-01   320.321
940       2025-05-01   320.580
941       2025-06-01   321.500
942       2025-07-01   322.132

[943 rows x 2 columns]

cpi['observation_date'] = pd.to_datetime(cpi['observation_date'])
cpi.set_index('observation_date', inplace = True)
print(cpi)

                  CPIAUCSL
observation_date          
1947-01-01          21.480
1947-02-01          21.620
1947-03-01          22.000
1947-04-01          22.000
1947-05-01          21.950
...                    ...
2025-03-01         319.615
2025-04-01         320.321
2025-05-01         320.580
2025-06-01         321.500
2025-07-01         322.132

[943 rows x 1 columns]

plt.figure(figsize=(12,5))
plt.plot(cpi.index, cpi['CPIAUCSL'], label='CPI')
plt.xlabel('Year')
plt.ylabel('CPI')
plt.title('Consumer Price Index (CPI)')
plt.show()

First, consider fitting a linear trend to the CPI data. This is given by

y_i = \beta_0 + \beta_1 x_i + \epsilon_i

where $y_i$ denotes the CPI for time index $i$ , and $x_i = i$ . The interpretation of $\beta_1$ is that it represents the numerical change in CPI for one month to the next. This is not a good model for CPI, because the relationship between CPI and data does not appear to be linear. Also CPI increases are usually reported in percentages (e.g., inflation is now at 2.8% etc.), and not in raw numerical increase. A better linear regression model here would use $y_i$ as the logarithm of CPI for time index $i$ :

y_i = \log(\text{CPI}_i) = \beta_0 + \beta_1 x_i + \epsilon_i

where $x_i = i$ as before. Now the interpretation of $\beta_1$ is that it represents the numerical change in $\log \text{CPI}_i$ for one month to the next. Note that

\log \text{CPI}_i - \log \text{CPI}_{i-1} = \log \frac{\text{CPI}_i}{\text{CPI}_{i-1}} \approx \frac{\text{CPI}_i}{\text{CPI}_{i-1}} - 1 = \frac{\text{CPI}_i - \text{CPI}_{i-1}}{\text{CPI}_{i-1}}

In other words, $100 \times \beta_1$ represents the percent change in CPI from one month to the next. If we want the percent change for one year (as opposed to one month), we can look at $12 \times 100 \times \beta_1$ . This number can be taken to represent an estimate of the historical inflation rate.

plt.figure(figsize=(12,4))
plt.plot(cpi.index, np.log(cpi['CPIAUCSL']), label='Logarithm of CPI')
plt.xlabel('Year')
plt.ylabel('Log(CPI)')
plt.title('Logarithm of Consumer Price Index (CPI)')
plt.show()

This model can be fit using the code below.

yvec = np.log(np.array(cpi['CPIAUCSL']))
n = cpi.shape[0]
xvec = np.arange(1, n+1)
X = sm.add_constant(xvec)
linreg = sm.OLS(yvec, X).fit()
print(linreg.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.966
Model:                            OLS   Adj. R-squared:                  0.966
Method:                 Least Squares   F-statistic:                 2.695e+04
Date:                Tue, 20 Oct 2026   Prob (F-statistic):               0.00
Time:                        19:55:28   Log-Likelihood:                 385.22
No. Observations:                 943   AIC:                            -766.4
Df Residuals:                     941   BIC:                            -756.7
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.9859      0.010    284.536      0.000       2.965       3.006
x1             0.0032   1.93e-05    164.171      0.000       0.003       0.003
==============================================================================
Omnibus:                     8570.174   Durbin-Watson:                   0.000
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               72.093
Skew:                          -0.026   Prob(JB):                     2.21e-16
Kurtosis:                       1.646   Cond. No.                     1.09e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.09e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

From this table, we see that the estimate of $\beta_0$ is 2.9859 and the estimate of $\beta_1$ is 0.0032. As per our interpretation, the month-to-month inflation rate is $100 \times \beta_1$ which is $0.32 \%$ . If we want the annual inflation rate, we would have to multiply this by 12 which gives:

historical_inflation_rate = 12 * 100 * linreg.params[1]
print(f'Historical annual inflation rate: {historical_inflation_rate:.3f}%')

Historical annual inflation rate: 3.794%

The regression output also gives an indication of the uncertainty associated with this inflation estimate:

conf_intervals_linreg = linreg.conf_int(alpha=0.05)
print(conf_intervals_linreg)
print(f'95% confidence interval for historical annual inflation rate (Linear Regression model): [{12*100*conf_intervals_linreg[1,0]:.3f}%, {   12*100*conf_intervals_linreg[1,1]:.3f}%]')

[[2.96527443 3.00646247]
 [0.003124   0.00319959]]
95% confidence interval for historical annual inflation rate (Linear Regression model): [3.749%, 3.840%]

This interval [3.749%, 3.84%] is extremely narrow for the historical inflation rate.

Below we plot the observed data along with the fitted values.

plt.figure(figsize=(12,4))
plt.plot(cpi.index, cpi['CPIAUCSL'], label='CPI')
plt.xlabel('Year')
plt.plot(cpi.index, np.exp(linreg.fittedvalues))
plt.show()

From this plot, it is clear that this regression model may not do a good job in predicting the values of CPI for the next few months. Below, we plot data vs fitted values on the logarithmic scale.

plt.figure(figsize=(12,4))
plt.plot(cpi.index, np.log(cpi['CPIAUCSL']), label='Logarithm of CPI')
plt.xlabel('Year')
plt.plot(cpi.index, linreg.fittedvalues)
plt.show()

Now, we shall formulate an alternative model resulting from placing the noise term in the evolution of the linear trend function, similar to how the Yule model was obtained from the sinusoid plus noise model.

The linear trend model can be written as:

\begin{align*} \log \text{CPI}_t = \ell(t) + \epsilon_t \end{align*}

(6)

where $\ell(t) = \beta_0 + \beta_1 t$ . Clearly $\ell(t) - \ell(t-1) = \beta_1$ (this is the discrete analogue of $\ell'(t) = \beta_1$ ). So an alternative linear plus noise model can be obtained by writing:

\begin{align*} \log \text{CPI}_t - \log \text{CPI}_{t-1} = \beta_1 + \epsilon_t. \end{align*}

(7)

This model is the same as:

\begin{align*} \log \text{CPI}_t - \log \text{CPI}_{t-1} \overset{\text{i.i.d}}{\sim} N(\beta_1, \sigma^2). \end{align*}

(8)

Thus, in this model, the estimate of $\beta_1$ is simply the mean of the values of $\log \text{CPI}_t - \log \text{CPI}_{t-1}$ . Below we compute the estimate of $\beta_1$ as well as the usual 95% confidence interval, compare this with the estimates and intervals for the historical inflation rate obtained via the previous model, and then argue that this is a better model than the previous one.

from scipy import stats

# Mean monthly log change and its usual 95% t interval.
dy = np.diff(yvec)
beta1_hat = np.mean(dy)
ci = stats.t.interval(0.95, len(dy) - 1, loc=beta1_hat, scale=stats.sem(dy))

comparison = pd.DataFrame({
    'Estimate (%)': [historical_inflation_rate, 1200 * beta1_hat],
    '95% lower (%)': [1200 * conf_intervals_linreg[1, 0], 1200 * ci[0]],
    '95% upper (%)': [1200 * conf_intervals_linreg[1, 1], 1200 * ci[1]]
}, index=['Log CPI vs time', 'Monthly log changes'])
comparison.round(3)

The estimate of $\beta_1$ here is somewhat smaller than the estimate of $\beta_1$ obtained from fitting linear regression (of $\log CPI_t$ on $t$ ). The confidence interval is also much more wider (about six times).

Below we plot the residuals in both the above models.

fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)

axes[0].plot(cpi.index, linreg.resid, color='black')
axes[0].set_title('Residuals: Log CPI vs time')

axes[1].plot(cpi.index[1:], dy - beta1_hat, color='blue')
axes[1].set_title('Residuals: Monthly log changes')

for ax in axes:
    ax.axhline(0, color='gray', linestyle='--')
    ax.set_ylabel('Residual')

axes[1].set_xlabel('Year')
plt.tight_layout()
plt.show()

From the first plot above, it is clear that $\log CPI_t = \ell(t) + \epsilon_t$ is not a good model for this dataset. In the second residual plot, there is some structure in the residuals but this suggests $\log CPI_t - \log CPI_{t-1} = \beta_1 + \epsilon_t$ is a much better model than the first one. Further improvements can be obtained by using AR models for $\log CPI_t - \log CPI_{t-1}$ as we shall see later.