How to obtain prediction intervals with statsmodels timeseries models?

Question

Is there a statsmodels API to retrieve prediction intervals from statsmodels timeseries models?

Currently, I'm manually calculating the prediction intervals using:

Here's my code. First, get some sample data ...

! python -c 'import datapackage' || pip install datapackage

%matplotlib inline

import datapackage

from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.api import SimpleExpSmoothing
import statsmodels.api as sm
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def get_data():
    # data licensed for non-commercial use only - https://datahub.io/core/bond-yields-uk-10y
    data_url = 'https://datahub.io/core/bond-yields-uk-10y/datapackage.json'

    resources = datapackage.Package(data_url).resources

    quarterly_csv_url = [pkg for pkg in resources if pkg.name == 'quarterly_csv'][0].descriptor['path']
    data = pd.read_csv(quarterly_csv_url)
    data = data.set_index('Date', drop=True).asfreq('Q')
    return data

Next, create a forecast and calculate the intervals:

data = get_data()
data = data[ data.index > '2005/']

fit = SimpleExpSmoothing(data).fit()
fcast = fit.forecast(1).rename('Forecast')
xhat = fcast.get_values()[0]

z = 1.96
sse = fit.sse
predint_xminus = xhat - z * np.sqrt(sse/len(data))
predint_xplus  = xhat + z * np.sqrt(sse/len(data))

Plot the intervals ...

plt.rcParams["figure.figsize"] = (20,5)

ax = data.plot(legend=True, title='British Goverment Bonds - 10y')
ax.set_xlabel('yield')

#
# 1-Step Prediction
#
prediction = pd.DataFrame( 
    data  = [ data.values[-1][0],  xhat ], 
    index = [ data.index[-1],      data.index[-1] + 1 ],
    columns = ['1-Step Predicted Rate']
)
_ = prediction.plot(ax=ax, color='black')

#
# upper 95% prediction interval
#
upper_pi_data = pd.DataFrame( 
    data  = [ xhat,           predint_xplus ], 
    index = [ data.index[-1], data.index[-1] + 1 ]
)
_ = upper_pi_data.plot(ax=ax, color='green', legend=False) 

#
# lower 95% prediction interval
#
lower_pi_data = pd.DataFrame( 
    data  = [ xhat,           predint_xminus ], 
    index = [ data.index[-1], data.index[-1] + 1 ]
)
_ = lower_pi_data.plot(ax=ax, color='green', legend=False) 

I've found similar questions, but not for timeseries models:

• confidence and prediction intervals with StatsModels

Answer 1

As long as you check the assumption that the residuals are uncorrelated and you don't go farther than one step ahead I think your prediction interval is valid. Note: I would use the standard deviation of the residuals. See section 3.5 in Forecasting Principles and Practice.

I'm pretty sure we need to place the model we are using into state space form according to Forecasting Principles and Practice for multi-step prediction intervals. See chapter 7.5 on exponential smoothing. The State Space Modeling for Local Linear Trend in statsmodels provides a working example. It doesn't look like there is anything out of the box to produce these intervals in statsmodels. I personally decided to use R to get my prediction intervals since the forecasting package provides these without a lot of additional effort.

Update: see comment below. Statsmodels now has state space representation for some exponential smoothing models.