<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 В файле записаны данные о заработной плате за каждый месяц с января 1993 по август 2016. Если хотите, можете дописать в конец ряда данные за следующие месяцы, если они уже опубликованы; найти эти данные можно, например, здесь.
Необходимо проанализировать данные, подобрать для них оптимальную прогнозирующую модель в классе ARIMA и построить прогноз на каждый месяц на два года вперёд от конца данных.
Review criteriaменьше 
Придерживайтесь стандартного алгоритма построения прогноза
Step1:  В рамках первичной визуалиции можно сразу отметить восходящий общий тренд. Сезонность с пиками в декабре и падением в январе(годовые премии). Рост дисперсии. Со всем этим нужно будет отдельно проанализировать
Step2:  Стабилизация дисперсии
Сделаем преобразование Бокса-Кокса для стабилизации дисперсии
Step3:  Стационарность
Попробуем сезонное дифференцирование - чтобы избавиться от тренда
Step4:  Это победа. Получили стационарный ряд. Lets rock
Step5:  SARIMA - наше всё
Посмотрим на ACF и PACF полученного ряда
Step6:  Начальные приближения
Step7:  Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Посмотрим, насколько хорошо модель описывает данные
Step8:  Модель описывает реальные данные очень даже хорошо
Прогноз
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn
import warnings
from itertools import product
import numpy as np 
def invboxcox(y,lmbda):
   if lmbda == 0:
      return(np.exp(y))
   else:
      return(np.exp(np.log(lmbda*y+1)/lmbda))
%matplotlib inline
salary = pd.read_csv('WAG_C_M.csv', sep=';', index_col=['month'], parse_dates=['month'], dayfirst=True)
# загрузили данные
salary.head()
salary.rename(columns={'WAG_C_M': 'salary_rub'}, inplace=True)
salary.salary_rub.plot(figsize=(15, 10), title='Average salary in Russia', fontsize=12);
plt.xlabel('month', fontsize=12)
plt.ylabel('average salary', fontsize=12)
plt.show()
Explanation: В файле записаны данные о заработной плате за каждый месяц с января 1993 по август 2016. Если хотите, можете дописать в конец ряда данные за следующие месяцы, если они уже опубликованы; найти эти данные можно, например, здесь.
Необходимо проанализировать данные, подобрать для них оптимальную прогнозирующую модель в классе ARIMA и построить прогноз на каждый месяц на два года вперёд от конца данных.
Review criteriaменьше 
Придерживайтесь стандартного алгоритма построения прогноза:
Визуальный анализ ряда
Стабилизация дисперсии (при необходимости)
Выбор порядка дифференцирования
Выбор начальных приближений для p,q,P,Q
Обучение и сравнение моделей-кандидатов, выбор победителя
Анализ остатков построенной модели, при необходимости — её модификация: если остатки получаются смещёнными,  прогноз нужно скорректировать на константу; если остатки нестационарны или автокоррелированы, можно попробовать расширить область, в которой подбираются значения p,q,P,Q.
Прогнозирование с помощью финальной модели.
End of explanation
# Проверка стационарности и STL-декомпозиция ряда:
plt.rcParams["figure.figsize"] = (10,15)
sm.tsa.seasonal_decompose(salary.salary_rub).plot( )
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(salary.salary_rub)[1])
# Гипотеза о стационарности критерием Дики- Фуллера не потверждается. Но подождите тут ведь еще тренд и сезонность
Explanation: В рамках первичной визуалиции можно сразу отметить восходящий общий тренд. Сезонность с пиками в декабре и падением в январе(годовые премии). Рост дисперсии. Со всем этим нужно будет отдельно проанализировать
End of explanation
plt.rcParams["figure.figsize"] = (10,8)
salary['salary_box'], lmbda = stats.boxcox(salary.salary_rub)
salary.salary_box.plot()
plt.ylabel('Transformed wine sales')
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(salary.salary_box)[1])
Explanation: Стабилизация дисперсии
Сделаем преобразование Бокса-Кокса для стабилизации дисперсии:
End of explanation
salary['salary_box_diff'] = salary.salary_box - salary.salary_box.shift(12)
sm.tsa.seasonal_decompose(salary.salary_box_diff[12:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(salary.salary_box_diff[12:])[1])
# Годовое диффернцирование не помогло, нужно еще раз
salary['salary_box_diff12_1'] = salary.salary_box_diff - salary.salary_box_diff.shift(1)
sm.tsa.seasonal_decompose(salary.salary_box_diff12_1[13:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(salary.salary_box_diff12_1[13:])[1])
Explanation: Стационарность
Попробуем сезонное дифференцирование - чтобы избавиться от тренда
End of explanation
salary.salary_box_diff = salary.salary_box_diff12_1
salary.drop('salary_box_diff12_1', axis=1, inplace=True)
Explanation: Это победа. Получили стационарный ряд. Lets rock
End of explanation
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(salary.salary_box_diff[13:].values.squeeze(), lags=48, ax=ax)
plt.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(salary.salary_box_diff[13:].values.squeeze(), lags=48, ax=ax)
plt.show()
Explanation: SARIMA - наше всё
Посмотрим на ACF и PACF полученного ряда
End of explanation
ps = range(0, 19)
d = 1
qs = range(0, 2)
Ps = range(0, 2)
D = 1
Qs = range(0, 1)
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(salary.salary_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())
# Лучшая модель:
print(best_model.summary())
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])
# Посмотрим на остатки модели:
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: Начальные приближения: Q = 0, q = 1, P = 1, p = 18
сезонных лагов со значимой корреляцией нет, значит, начальное приближение Q = 0
End of explanation
salary['model'] = invboxcox(best_model.fittedvalues, lmbda)
plt.rcParams["figure.figsize"] = (10,8)
salary.salary_rub.plot()
salary.model[13:].plot(color='r')
plt.ylabel('Salary in Russia')
plt.show()
Explanation: Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Посмотрим, насколько хорошо модель описывает данные:
End of explanation
from datetime import datetime
import datetime
from dateutil.relativedelta import *
salary2 = salary[['salary_rub']]
date_list = [datetime.datetime.strptime("2016-09-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,36)]
future = pd.DataFrame(index=date_list, columns=salary2.columns)
salary2 = pd.concat([salary2, future])
salary2['forecast'] = invboxcox(best_model.predict(start=284, end=325), lmbda)
salary2.salary_rub.plot()
salary2.forecast.plot(color='r')
plt.ylabel('Average salary in Russia (rub)')
plt.show()
Explanation: Модель описывает реальные данные очень даже хорошо
Прогноз
End of explanation
<END_TASK>