from __future__ import print_function
%matplotlib inline
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.statespace import structural
import matplotlib.pyplot as plt
# Get dataset
nile = sm.datasets.nile.data.load_pandas().data
nile.index = pd.date_range('1871', '1970', freq='AS')
# Create exogenous array for intervention in 1899
# and for pulse in 1913
nile['intervention'] = 0
nile.ix['1899':, 'intervention'] = 1
nile['pulse'] = 0
nile.ix['1913', 'pulse'] = 1
nile['volume'].plot(figsize=(10,3));
mod_ll = structural.UnobservedComponents(
nile['volume'], exog=nile[['intervention', 'pulse']],
level=True, stochastic_level=True,
)
res_ll = mod_ll.fit()
res_ll = mod_ll.smooth()
print(res_ll.summary())
Statespace Model Results ================================================================================ Dep. Variable: volume No. Observations: 100 Model: UnobservedComponents Log Likelihood -617.215 Date: Wed, 03 Jun 2015 AIC 1242.431 Time: 10:59:51 BIC 1252.852 Sample: 01-01-1871 HQIC 1246.648 - 01-01-1970 ==================================================================================== coef std err z P>|z| [95.0% Conf. Int.] ------------------------------------------------------------------------------------ sigma2.irregular 1.455e+04 2064.796 7.045 0.000 1.05e+04 1.86e+04 sigma2.level 8.351e-05 nan nan nan nan nan beta.1 -241.5867 23.438 -10.308 0.000 -287.524 -195.650 beta.2 -399.5548 121.299 -3.294 0.001 -637.296 -161.814 ====================================================================================
fig, axes = plt.subplots(2, 1, figsize=(10,6))
dates = nile.index._mpl_repr()
# Predictions, confidence intervals
ax = axes[0]
ax.plot(dates, nile['volume'], 'k.')
ax.plot(dates, res_ll.forecasts[0])
ax.set_title('Filtered level');
# Filtered level
ax = axes[1]
ax.plot(dates, nile['volume'], 'k.')
ax.plot(dates, res_ll.smoothed_state[0] + res_ll.obs_intercept[0])
ax.set_title('Smoothed level');
# Create the model
mod_lc = structural.UnobservedComponents(
nile['volume'], exog=nile[['intervention','pulse']],
irregular=False,
level=True, stochastic_level=False,
cycle=True, stochastic_cycle=True, damped_cycle=True
)
# Fit the model
res_lc = mod_lc.fit()
res_lc = mod_lc.smooth()
print(res_lc.summary())
Statespace Model Results ================================================================================ Dep. Variable: volume No. Observations: 100 Model: UnobservedComponents Log Likelihood -604.025 Date: Wed, 03 Jun 2015 AIC 1218.050 Time: 10:59:51 BIC 1231.076 Sample: 01-01-1871 HQIC 1223.322 - 01-01-1970 =================================================================================== coef std err z P>|z| [95.0% Conf. Int.] ----------------------------------------------------------------------------------- sigma2.cycle 1.435e+04 2062.830 6.955 0.000 1.03e+04 1.84e+04 frequency.cycle 0.5236 1.270 0.412 0.680 -1.966 3.013 damping.cycle 0.1370 0.112 1.228 0.219 -0.082 0.356 beta.1 -241.1206 30.793 -7.830 0.000 -301.473 -180.768 beta.2 -380.6163 121.617 -3.130 0.002 -618.981 -142.251 ===================================================================================
fig, axes = plt.subplots(2, 1, figsize=(10,6))
dates = nile.index._mpl_repr()
# Plot the level
ax = axes[0]
ax.plot(dates, nile['volume'], 'k.')
ax.plot(dates, res_lc.filtered_state[0] + res_lc.obs_intercept[0], 'r', label='Filtered')
ax.plot(dates, res_lc.smoothed_state[0] + res_lc.obs_intercept[0], 'k', label='Smoothed')
ax.legend()
ax.set_title('Level');
# Plot the cycle
ax = axes[1]
ax.plot(dates, res_lc.filtered_state[1], 'r', label='Filtered')
ax.plot(dates, res_lc.smoothed_state[1], 'k', label='Smoothed')
ax.legend()
ax.set_title('Cycle');
fig, axes = plt.subplots(2, 1, figsize=(10,9))
dates = nile.index._mpl_repr()
# Plot the level
ax = axes[0]
ax.plot(dates, nile['volume'], 'k.')
ax.plot(dates, res_lc.filtered_state[0] + res_lc.obs_intercept[0])
ax.set_title('Smoothed level');
# Plot the cycle
ax = axes[1]
ax.plot(dates, res_lc.filtered_state[1])
ax.set_title('Filtered cycle');