title: How to obtain prediction intervals with statsmodels timeseries models?
body: <p>Is there a statsmodels API to retrieve prediction intervals from statsmodels timeseries models?</p>
  
  <p>Currently, I'm manually calculating the prediction intervals using:</p>
  
  <p><a href="https://i.stack.imgur.com/wdwMH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/wdwMH.png" alt="enter image description here"></a></p>
  
  <p>Here's my code.  First, get some sample data ...</p>
  
  <pre class="lang-py prettyprint-override"><code>! 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
  </code></pre>
  
  <p>Next, create a forecast and calculate the intervals:</p>
  
  <pre class="lang-py prettyprint-override"><code>data = get_data()
  data = data[ data.index &gt; '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))
  </code></pre>
  
  <p>Plot the intervals ...</p>
  
  <pre class="lang-py prettyprint-override"><code>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) 
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/Omx7g.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Omx7g.png" alt="enter image description here"></a></p>
  
  <p>I've found similar questions, but not for timeseries models:</p>
  
  <ul>
  <li><a href="https://stackoverflow.com/questions/17559408/confidence-and-prediction-intervals-with-statsmodels">confidence and prediction intervals with StatsModels</a></li>
  </ul>