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()
print(res_ll.summary())
res_ll.plot_diagnostics(figsize=(11,6));
Statespace Model Results ================================================================================ Dep. Variable: volume No. Observations: 100 Model: UnobservedComponents Log Likelihood -617.221 Date: Tue, 02 Jun 2015 AIC 1242.441 Time: 16:00:02 BIC 1252.862 Sample: 01-01-1871 HQIC 1246.659 - 01-01-1970 ==================================================================================== coef std err z P>|z| [95.0% Conf. Int.] ------------------------------------------------------------------------------------ sigma2.irregular 1.455e+04 2067.494 7.040 0.000 1.05e+04 1.86e+04 sigma2.level 1.988e-05 nan nan nan nan nan beta.1 -242.4507 22.796 -10.635 0.000 -287.131 -197.771 beta.2 -387.6257 121.323 -3.195 0.001 -625.415 -149.836 =================================================================================== Ljung-Box (Q): 41.48 Jarque-Bera (JB): 0.82 Prob(Q): 0.41 Prob(JB): 0.66 Heteroskedasticity (H): 0.88 Skew: 0.18 Prob(H) (two-sided): 0.72 Kurtosis: 2.74 ===================================================================================
fig, axes = plt.subplots(2, 1, figsize=(10,6))
dates = nile.index._mpl_repr()
# Filtered level
ax = axes[0]
ax.plot(dates, nile['volume'], 'k.')
ax.plot(dates, res_ll.filtered_state[0] + res_ll.obs_intercept[0])
ax.set_title('Filtered level');
# Residuals
ax = axes[1]
ax.plot(dates, nile['volume'] - res_ll.filtered_state[0]- res_lc.obs_intercept[0])
ax.set_title('Filtered residuals');
# 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()
print(res_lc.summary())
res_lc.plot_diagnostics(figsize=(11,6));
Statespace Model Results ================================================================================ Dep. Variable: volume No. Observations: 100 Model: UnobservedComponents Log Likelihood -604.025 Date: Tue, 02 Jun 2015 AIC 1218.050 Time: 15:53:11 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 =================================================================================== Ljung-Box (Q): 37.78 Jarque-Bera (JB): 0.77 Prob(Q): 0.57 Prob(JB): 0.68 Heteroskedasticity (H): 0.80 Skew: 0.18 Prob(H) (two-sided): 0.52 Kurtosis: 2.76 ===================================================================================
fig, axes = plt.subplots(3, 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('Filtered level');
# Plot the cycle
ax = axes[1]
ax.plot(dates, res_lc.filtered_state[1])
ax.set_title('Filtered cycle');
# Plot the filtered residuals
ax = axes[2]
ax.plot(dates, nile['volume'] - res_lc.filtered_state[0] - res_lc.filtered_state[1] - res_lc.obs_intercept[0])
ax.set_title('Filtered residuals')
ax.set_ylim(-300,400);