状態空間法による時系列解析 statespace¶

# Load the statsmodels api
import statsmodels.api as sm

# Load your dataset
endog = pd.read_csv('your/dataset/here.csv')

# We could fit an AR(2) model, described above
mod_ar2 = sm.tsa.SARIMAX(endog, order=(2,0,0))
# Note that mod_ar2 is an instance of the SARIMAX class

# Fit the model via maximum likelihood
res_ar2 = mod_ar2.fit()
# Note that res_ar2 is an instance of the SARIMAXResults class

# Show the summary of results
print(res_ar2.summary())

# We could also fit a more complicated model with seasonal components.
# As an example, here is an SARIMA(1,1,1) x (0,1,1,4):
mod_sarimax = sm.tsa.SARIMAX(endog, order=(1,1,1),
                             seasonal_order=(0,1,1,4))
res_sarimax = mod_sarimax.fit()

# Show the summary of results
print(res_sarimax.summary())

# Load the statsmodels api
import statsmodels.api as sm

# Load your dataset
endog = pd.read_csv('your/dataset/here.csv')

# Fit a local level model
mod_ll = sm.tsa.UnobservedComponents(endog, 'local level')
# Note that mod_ll is an instance of the UnobservedComponents class

# Fit the model via maximum likelihood
res_ll = mod_ll.fit()
# Note that res_ll is an instance of the UnobservedComponentsResults class

# Show the summary of results
print(res_ll.summary())

# Show a plot of the estimated level and trend component series
fig_ll = res_ll.plot_components()

# We could further add a damped stochastic cycle as follows
mod_cycle = sm.tsa.UnobservedComponents(endog, 'local level', cycle=True,
                                        damped_cycle=True,
                                        stochastic_cycle=True)
res_cycle = mod_cycle.fit()

# Show the summary of results
print(res_cycle.summary())

# Show a plot of the estimated level, trend, and cycle component series
fig_cycle = res_cycle.plot_components()

# Load the statsmodels api
import statsmodels.api as sm

# Load your (multivariate) dataset
endog = pd.read_csv('your/dataset/here.csv')

# Fit a local level model
mod_var1 = sm.tsa.VARMAX(endog, order=(1,0))
# Note that mod_var1 is an instance of the VARMAX class

# Fit the model via maximum likelihood
res_var1 = mod_var1.fit()
# Note that res_var1 is an instance of the VARMAXResults class

# Show the summary of results
print(res_var1.summary())

# Construct impulse responses
irfs = res_ll.impulse_responses(steps=10)

# Load the statsmodels api
import statsmodels.api as sm

# Load your dataset
endog = pd.read_csv('your/dataset/here.csv')

# Create a dynamic factor model
mod_dfm = sm.tsa.DynamicFactorMQ(endog, k_factors=1, factor_order=2)
# Note that mod_dfm is an instance of the DynamicFactorMQ class

# Fit the model via maximum likelihood, using the EM algorithm
res_dfm = mod_dfm.fit()
# Note that res_dfm is an instance of the DynamicFactorMQResults class

# Show the summary of results
print(res_ll.summary())

# Show a plot of the r^2 values from regressions of
# individual estimated factors on endogenous variables.
fig_dfm = res_ll.plot_coefficients_of_determination()

# Load the statsmodels api
import statsmodels.api as sm

# Load your dataset
endog = pd.read_csv('your/dataset/here.csv')

# Simple exponential smoothing, denoted (A,N,N)
mod_ses = sm.tsa.statespace.ExponentialSmoothing(endog)
res_ses = mod_ses.fit()

# Holt's linear method, denoted (A,A,N)
mod_h = sm.tsa.statespace.ExponentialSmoothing(endog, trend=True)
res_h = mod_h.fit()

# Damped trend model, denoted (A,Ad,N)
mod_dt = sm.tsa.statespace.ExponentialSmoothing(endog, trend=True,
                                                damped_trend=True)
res_dt = mod_dt.fit()

# Holt-Winters' trend and seasonality method, denoted (A,A,A)
# (assuming that `endog` has a seasonal periodicity of 4, for example if it
# is quarterly data).
mod_hw = sm.tsa.statespace.ExponentialSmoothing(endog, trend=True,
                                                seasonal=4)
res_hw = mod_hw.fit()

import numpy as np
from scipy.signal import lfilter
import statsmodels.api as sm

# True model parameters
nobs = int(1e3)
true_phi = np.r_[0.5, -0.2]
true_sigma = 1**0.5

# Simulate a time series
np.random.seed(1234)
disturbances = np.random.normal(0, true_sigma, size=(nobs,))
endog = lfilter([1], np.r_[1, -true_phi], disturbances)

# Construct the model
class AR2(sm.tsa.statespace.MLEModel):
    def __init__(self, endog):
        # Initialize the state space model
        super(AR2, self).__init__(endog, k_states=2, k_posdef=1,
                                  initialization='stationary')

        # Setup the fixed components of the state space representation
        self['design'] = [1, 0]
        self['transition'] = [[0, 0],
                                  [1, 0]]
        self['selection', 0, 0] = 1

    # Describe how parameters enter the model
    def update(self, params, transformed=True, **kwargs):
        params = super(AR2, self).update(params, transformed, **kwargs)

        self['transition', 0, :] = params[:2]
        self['state_cov', 0, 0] = params[2]

    # Specify start parameters and parameter names
    @property
    def start_params(self):
        return [0,0,1]  # these are very simple

# Create and fit the model
mod = AR2(endog)
res = mod.fit()
print(res.summary())

# Load in the example macroeconomic dataset
dta = sm.datasets.macrodata.load_pandas().data
# Make sure we have an index with an associated frequency, so that
# we can refer to time periods with date strings or timestamps
dta.index = pd.date_range('1959Q1', '2009Q3', freq='QS')

# Step 1: construct an SARIMAX model for US inflation data
model = sm.tsa.SARIMAX(dta.infl, order=(4, 0, 0), trend='c')

# Step 2: fit the model's parameters by maximum likelihood
results = model.fit()

# Step 3: explore / use results

# - Print a table summarizing estimation results
print(results.summary())

# - Print only the estimated parameters
print(results.params)

# - Create diagnostic figures based on standardized residuals:
#   (1) time series graph
#   (2) histogram
#   (3) Q-Q plot
#   (4) correlogram
results.plot_diagnostics()

# - Examine diagnostic hypothesis tests
# Jarque-Bera: [test_statistic, pvalue, skewness, kurtosis]
print(results.test_normality(method='jarquebera'))
# Goldfeld-Quandt type test: [test_statistic, pvalue]
print(results.test_heteroskedasticity(method='breakvar'))
# Ljung-Box test: [test_statistic, pvalue] for each lag
print(results.test_serial_correlation(method='ljungbox'))

# - Forecast the next 4 values
print(results.forecast(4))

# - Forecast until 2020Q4
print(results.forecast('2020Q4'))

# - Plot in-sample dynamic prediction starting in 2005Q1
#   and out-of-sample forecasts until 2010Q4 along with
#   90% confidence intervals
predict_results = results.get_prediction(start='2005Q1', end='2010Q4', dynamic=True)
predict_df = predict_results.summary_frame(alpha=0.10)
fig, ax = plt.subplots()
predict_df['mean'].plot(ax=ax)
ax.fill_between(predict_df.index, predict_df['mean_ci_lower'],
                predict_df['mean_ci_upper'], alpha=0.2)

# - Simulate two years of new data after the end of the sample
print(results.simulate(8, anchor='end'))

# - Impulse responses for two years
print(results.impulse_responses(8))

fig, axes = plt.subplots(3, 1, figsize=(8, 8))

# Retrieve monthly retail sales for clothing
from pandas_datareader.data import DataReader
clothing = DataReader('MRTSSM4481USN', 'fred', start='1992').asfreq('MS')['MRTSSM4481USN']

# Construct a local level + seasonal model
model = sm.tsa.UnobservedComponents(clothing, 'llevel', seasonal=12)
results = model.fit()

# Plot the data, the level, and seasonal
clothing.plot(ax=axes[0])
results.states.smoothed['level'].plot(ax=axes[1])
results.states.smoothed['seasonal'].plot(ax=axes[2])

# Load in the example macroeconomic dataset
dta = sm.datasets.macrodata.load_pandas().data
# Make sure we have an index with an associated frequency, so that
# we can refer to time periods with date strings or timestamps
dta.index = pd.date_range('1959Q1', '2009Q3', freq='QS')

# Separate inflation data into a training and test dataset
training_endog = dta['infl'].iloc[:-1]
test_endog = dta['infl'].iloc[-1:]

# Fit an SARIMAX model for inflation
training_model = sm.tsa.SARIMAX(training_endog, order=(4, 0, 0))
training_results = training_model.fit()

# Extend the results to the test observations
test_results = training_results.extend(test_endog)

# Print the sum of squared errors in the test sample,
# based on parameters computed using only the training sample
print(test_results.sse)

# Construct a model
model = sm.tsa.SARIMAX(endog, order=(1, 0, 0))

# To find out the parameter names, use:
print(model.param_names)

# Fit the model with a fixed value for the AR(1) coefficient:
results = model.fit_constrained({'ar.L1': 0.5})

# Construct a model
model = sm.tsa.SARIMAX(endog, order=(1, 0, 0))

# Fit the model with a fixed value for the AR(1) coefficient using the
# context manager
with model.fix_params({'ar.L1': 0.5}):
    results = model.fit()