How to use the statsmodels.api.tsa function in statsmodels

maparams = np.array([.65, .35])
# The conventions of the arma_generate function require that we specify a
# 1 for the zero-lag of the AR and MA parameters and that the AR parameters
# be negated.
arparams = np.r_[1, -arparams]
maparam = np.r_[1, maparams]
nobs = 250
y = arma_generate_sample(arparams, maparams, nobs)
plt.figure()
plt.plot(y)

# Now, optionally, we can add some dates information. For this example,
# we'll use a pandas time series.
dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)
y = Series(y, index=dates)
arma_mod = sm.tsa.ARMA(y, order=(2, 2), freq='M')
arma_res = arma_mod.fit(trend='nc', disp=-1)
print(arma_res.params)

plt.show()

# unrestricted_model = {
#     'irregular': True, 'level': True, 'stochastic_level': False,
# 'trend': True, 'stochastic_trend': True,
#     'cycle': True, 'damped_cycle': True, 'stochastic_cycle': True
# }

# We now fit the following models:
#
# 1. Output, unrestricted model
# 2. Prices, unrestricted model
# 3. Prices, restricted model
# 4. Money, unrestricted model
# 5. Money, restricted model

# Output
output_mod = sm.tsa.UnobservedComponents(dta['US GNP'], **unrestricted_model)
output_res = output_mod.fit(method='powell', disp=False)

# Prices
prices_mod = sm.tsa.UnobservedComponents(dta['US Prices'],
                                         **unrestricted_model)
prices_res = prices_mod.fit(method='powell', disp=False)

prices_restricted_mod = sm.tsa.UnobservedComponents(dta['US Prices'],
                                                    **restricted_model)
prices_restricted_res = prices_restricted_mod.fit(method='powell', disp=False)

# Money
money_mod = sm.tsa.UnobservedComponents(dta['US monetary base'],
                                        **unrestricted_model)
money_res = money_mod.fit(method='powell', disp=False)

cpi = load_pandas().data['cpi']
dates = dates_from_range('1959q1', '2009q3')
cpi.index = dates

res = ARIMA(cpi, (1, 1, 1), freq='Q').fit()
print(res.summary())

# we can look at the series
cpi.diff().plot()

# maybe logs are better
log_cpi = np.log(cpi)

# check the ACF and PCF plots
acf, confint_acf = sm.tsa.acf(log_cpi.diff().values[1:], confint=95)
# center the confidence intervals about zero
#confint_acf -= confint_acf.mean(1)[:, None]
pacf = sm.tsa.pacf(log_cpi.diff().values[1:], method='ols')
# confidence interval is now an option to pacf
from scipy import stats
confint_pacf = stats.norm.ppf(1 - .025) * np.sqrt(1 / 202.)

fig = plt.figure()
ax = fig.add_subplot(121)
ax.set_title('Autocorrelation')
ax.plot(range(41), acf, 'bo', markersize=5)
ax.vlines(range(41), 0, acf)
ax.fill_between(range(41), confint_acf[:, 0], confint_acf[:, 1], alpha=.25)
fig.tight_layout()
ax = fig.add_subplot(122, sharey=ax)
ax.vlines(range(41), 0, pacf)

cpi.index = dates

res = ARIMA(cpi, (1, 1, 1), freq='Q').fit()
print(res.summary())

# we can look at the series
cpi.diff().plot()

# maybe logs are better
log_cpi = np.log(cpi)

# check the ACF and PCF plots
acf, confint_acf = sm.tsa.acf(log_cpi.diff().values[1:], confint=95)
# center the confidence intervals about zero
#confint_acf -= confint_acf.mean(1)[:, None]
pacf = sm.tsa.pacf(log_cpi.diff().values[1:], method='ols')
# confidence interval is now an option to pacf
from scipy import stats
confint_pacf = stats.norm.ppf(1 - .025) * np.sqrt(1 / 202.)

fig = plt.figure()
ax = fig.add_subplot(121)
ax.set_title('Autocorrelation')
ax.plot(range(41), acf, 'bo', markersize=5)
ax.vlines(range(41), 0, acf)
ax.fill_between(range(41), confint_acf[:, 0], confint_acf[:, 1], alpha=.25)
fig.tight_layout()
ax = fig.add_subplot(122, sharey=ax)
ax.vlines(range(41), 0, pacf)
ax.plot(range(41), pacf, 'bo', markersize=5)
ax.fill_between(range(41), -confint_pacf, confint_pacf, alpha=.25)

# Instantiate the model

ar_model = sm.tsa.AR(endog, freq='A')
pandas_ar_res = ar_model.fit(maxlag=9, method='mle', disp=-1)


# Out-of-sample prediction

pred = pandas_ar_res.predict(start='2005', end='2015')
print(pred)


# ## Using explicit dates

ar_model = sm.tsa.AR(data.endog, dates=dates, freq='A')
ar_res = ar_model.fit(maxlag=9, method='mle', disp=-1)
pred = ar_res.predict(start='2005', end='2015')
print(pred)


# This just returns a regular array, but since the model has date information attached, you can get the prediction dates in a roundabout way.

print(ar_res.data.predict_dates)

# The conventions of the arma_generate function require that we specify a 1 for the zero-lag of the AR and MA parameters and that the AR parameters be negated.

arparams = np.r_[1, -arparams]
maparam = np.r_[1, maparams]
nobs = 250
y = arma_generate_sample(arparams, maparams, nobs)


#  Now, optionally, we can add some dates information. For this example, we'll use a pandas time series.

import pandas as pd
dates = sm.tsa.datetools.dates_from_range('1980m1', length=nobs)
y = pd.Series(y, index=dates)
arma_mod = sm.tsa.ARMA(y, order=(2,2))
arma_res = arma_mod.fit(trend='nc', disp=-1)

It produces false positives on non-stationary series so Augmented
    Dickey-Fuller test applied to check for stationarity.
    """

    hour_ago = time() - 3600
    ten_minutes_ago = time() - 600
    reference = scipy.array([x[1] for x in timeseries if x[0] >= hour_ago and x[0] < ten_minutes_ago])
    probe = scipy.array([x[1] for x in timeseries if x[0] >= ten_minutes_ago])

    if reference.size < 20 or probe.size < 20:
        return False

    ks_d, ks_p_value = scipy.stats.ks_2samp(reference, probe)

    if ks_p_value < 0.05 and ks_d > 0.5:
        adf = sm.tsa.stattools.adfuller(reference, 10)
        if adf[1] < 0.05:
            return True

    return False

# folder with data
data_folder = '../../Data/Chapter07/'

# colors 
colors = ['#FF6600', '#000000', '#29407C', '#660000']

# read the data
riverFlows = pd.read_csv(data_folder + 'combined_flow.csv', 
    index_col=0, parse_dates=[0])

# autocorrelation function
acf = {}    # to store the results
f = {}

for col in riverFlows.columns:
    acf[col] = sm.tsa.stattools.acf(riverFlows[col])

# partial autocorrelation function
pacf = {}

for col in riverFlows.columns:
    pacf[col] = sm.tsa.stattools.pacf(riverFlows[col])

# periodogram (spectral density)
sd = {}

for col in riverFlows.columns:
    sd[col] = sm.tsa.stattools.periodogram(riverFlows[col])

# plot the data
fig, ax = plt.subplots(2, 3) # 2 rows and 3 columns

def RstarStage1(self):

        print('===== STAGE 1 =====')

        stage = 1

        # Data must start 4 quarters before the estimation period
        T = self.logGDP.shape[0] - 4

        # Original output gap estimate (linear time trend)
        x_og = np.concatenate((np.ones((T, 1)), np.arange(1, T + 1, 1).reshape((T, 1))), axis=1)
        y_og = self.logGDP[4:].values.reshape((T, 1))
        output_gap = (y_og - x_og.dot(np.linalg.solve(x_og.T.dot(x_og), x_og.T.dot(y_og)))) * 100

        # Initialization of state vector for the Kalman filter using HP trend of log output
        logGDP_HP_cycle, logGDP_HP_trend = sm.tsa.filters.hpfilter(self.logGDP, 1600)
        g_pot = np.flip(logGDP_HP_trend[0:3].values, axis=0).reshape((3, 1))
        xi_00 = 100*g_pot

        # IS curve
        y_is = output_gap[2: T]
        # y_is = logGDP_HP_cycle.values[6:].reshape((T - 2, 1))
        x_is = np.concatenate((output_gap[1:T-1], output_gap[0:T-2]), axis=1)
        b_is = np.linalg.solve(x_is.T.dot(x_is), x_is.T.dot(y_is))  # b stands for beta
        r_is = y_is - x_is.dot(b_is)  # r stands for residuals
        s_is = np.sqrt(r_is.T.dot(r_is) / (r_is.shape[0] - x_is.shape[1]))[0, 0]  # s stands for standard deviation of residuals (unbiased estimate)

        # Phillips Curve
        y_ph = self.inflation[4: T].values.reshape((T-4, 1))  # RESHAPE T-4
        x_ph = np.concatenate((self.inflation[3: T-1].values.reshape((T-4, 1)),
                               (self.inflation[2: T-2].values.reshape((T-4, 1)) + self.inflation[1: T-3].values.reshape((T-4, 1)) + self.inflation[0: T-4].values.reshape((T-4, 1))) / 3,
                               output_gap[3: T-1]), axis=1)

def arma_sample(n, ar, ma, verbose=False):
    '''
    Generate an arma model from lags by wrapping statsmodels.
    '''
    zero_lag_ar = np.r_[1, -np.array(ar)]
    zero_lag_ma = np.r_[1, np.array(ma)]

    arma_process = sm.tsa.ArmaProcess(zero_lag_ar, zero_lag_ma)
    sample = arma_process.generate_sample(n)
    if verbose:
        print('AR roots: ', arma_process.arroots)
        print('MA roots: ', arma_process.maroots)

    return sample