<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 свойства ряда
Step1:  тренд имеет простую, легко объяснимую инфляцией, структуру
сезонный профиль четко выражен
следующий шаг - стабилизация дисперсии
Step2:  И визуально, и воспользовашись критерием Дики-Фуллера, мы можем понять, что ряд по-прежнему не является стационарным. Но размах сезонных колебаний значительно уменьшился. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность
Step3:  Критерий Дики-Фуллера отвергает гипотезу нестационарности ряда на уровне значимости 0.05, но мы по-прежнему видим тренд. Попробуем добавить еще обычное дифференцирование
Step4:  Гипотеза нестационарности с уверенностью отвергается, и визуально ряд выглядит лучше. На трендовой компоненте не видно никакого систематического поведения - график колеблется в районе 0. 
Перейдем к посмотрению модели. Для этого посмотрим на ACF и PACF полученного ряда
Step5:  максимальный сезонный лаг отличающийся от нуля отсутствует => Q = 0
последний несезонный лаг, при котором автокорреляция значима, 1 => q = 1
последний сезонный лаг, при котором частичная автокорреляция значима, 12 => P = 1
номер последнего несезонного лага, при котором частичная автокорреляция значима, 1 => q = 1 (значения остальных значимых лагов превышают сезонность, 12)
Step6:  переберем все возможные значения параметров P, p, Q, q
Step7:  Лучшая модель (минимальный aic = 38.9, parameters = (1, 0, 1, 0))
Step8:  Остатки
Step9:  Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой).
Посмотрим, насколько хорошо модель описывает данные
Step10:  Прогноз
<ASSISTANT_TASK:>
Python Code:
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.WAG_C_M).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.WAG_C_M)[1])
Explanation: свойства ряда:
* повышающийся тренд
* годовая сезонность
* автокоррелированность
* в начале ряда размах сезонных колебаний значительно меньше, чем в конце => дисперсия нестационарна
Можно с уверенностью утверждать, что ряд не обладает свойствами стационарности, но тем не менее воспользуемся критерием Дики-Фуллера, а также проведем STL-декомпозицию ряда:
End of explanation
wages['wages_box'], lmbda = stats.boxcox(wages.WAG_C_M)
plt.figure(figsize(15,7))
wages.wages_box.plot()
plt.ylabel('Transformed wages')
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box)[1])
Explanation: тренд имеет простую, легко объяснимую инфляцией, структуру
сезонный профиль четко выражен
следующий шаг - стабилизация дисперсии:
End of explanation
wages['wages_box_diff'] = wages.wages_box - wages.wages_box.shift(12)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.wages_box_diff[12:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box_diff[12:])[1])
Explanation: И визуально, и воспользовашись критерием Дики-Фуллера, мы можем понять, что ряд по-прежнему не является стационарным. Но размах сезонных колебаний значительно уменьшился. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:
End of explanation
wages['wages_box_diff2'] = wages.wages_box_diff - wages.wages_box_diff.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(wages.wages_box_diff2[13:]).plot()   
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(wages.wages_box_diff2[13:])[1])
Explanation: Критерий Дики-Фуллера отвергает гипотезу нестационарности ряда на уровне значимости 0.05, но мы по-прежнему видим тренд. Попробуем добавить еще обычное дифференцирование:
End of explanation
plt.figure(figsize(15,8))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(wages.wages_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(wages.wages_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()
Explanation: Гипотеза нестационарности с уверенностью отвергается, и визуально ряд выглядит лучше. На трендовой компоненте не видно никакого систематического поведения - график колеблется в районе 0. 
Перейдем к посмотрению модели. Для этого посмотрим на ACF и PACF полученного ряда:
End of explanation
ps = range(0, 2)
d=1
qs = range(0, 2)
Ps = range(0, 2)
D=1
Qs = range(0, 1)
Explanation: максимальный сезонный лаг отличающийся от нуля отсутствует => Q = 0
последний несезонный лаг, при котором автокорреляция значима, 1 => q = 1
последний сезонный лаг, при котором частичная автокорреляция значима, 12 => P = 1
номер последнего несезонного лага, при котором частичная автокорреляция значима, 1 => q = 1 (значения остальных значимых лагов превышают сезонность, 12)
End of explanation
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
    #try except нужен, потому что на некоторых наборах параметров модель не обучается
    try:
        model=sm.tsa.statespace.SARIMAX(wages.wages_box, order=(param[0], d, param[1]), 
                                        seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
    #выводим параметры, на которых модель не обучается и переходим к следующему набору
    except ValueError:
        print('wrong parameters:', param)
        continue
    aic = model.aic
    #сохраняем лучшую модель, aic, параметры
    if aic < best_aic:
        best_model = model
        best_aic = aic
        best_param = param
    results.append([param, model.aic])
    
warnings.filterwarnings('default')
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
Explanation: переберем все возможные значения параметров P, p, Q, q:
End of explanation
print(best_model.summary())
Explanation: Лучшая модель (минимальный aic = 38.9, parameters = (1, 0, 1, 0)):
End of explanation
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid[13:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])
Explanation: Остатки:
End of explanation
def invboxcox(y,lmbda):
   if lmbda == 0:
      return(np.exp(y))
   else:
      return(np.exp(np.log(lmbda*y+1)/lmbda))
wages['model'] = invboxcox(best_model.fittedvalues, lmbda)
plt.figure(figsize(15,7))
wages.WAG_C_M.plot()
wages.model[13:].plot(color='r')
plt.ylabel('wages')
pylab.show()
Explanation: Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой).
Посмотрим, насколько хорошо модель описывает данные:
End of explanation
wages2 = wages[['WAG_C_M']]
date_list = [datetime.datetime.strptime("2017-07-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,24)]
future = pd.DataFrame(index=date_list, columns= wages2.columns)
wages2 = pd.concat([wages2, future])
wages2['forecast'] = invboxcox(best_model.predict(start=294, end=317), lmbda)
plt.figure(figsize(15,7))
wages2.WAG_C_M.plot()
wages2.forecast.plot(color='r')
plt.ylabel('wages')
pylab.show()
Explanation: Прогноз:
End of explanation
<END_TASK>