<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Because the p-value is larger than 0.05, the moving average is not constant over time and the null hypothesis of the Dickey-Fuller test cannot be rejected. This shows that the weekly time series is not stationary. Before you can apply ARIMA models for forecasting, you need to transform this time series into a stationary time series.
Step1:  Because the p-value is now smaller than 0.05, we have necessary evidence to reject the null hypothesis of the Dickey-Fuller test. This shows that the time series is stationary at log scale. We can apply ARIMA model forcasting on this log scaled data.
Remove trend and seasonality with decomposition¶
Step2:  Now that you have stationarized your time series, you could go on and model residuals (fit lines between values in the plot). However, as the patterns for the trend and seasonality information extracted from the series that are plotted after decomposition are still not consistent and cannot be scaled back to the original values, you cannot use this approach to create reliable forecasts.
In case your time series does show a strong and consistent seasonal trend, here is a good article that describes the use of the SARIMAX model in Python.
Remove trend and seasonality with differencing
Step3:  Find optimal parameters and build an ARIMA (auto regressive integrated moving average) model
To apply an ARIMA model to your time series, you need to find optimal values for the following three model parameters (p,d,q)
Step4:  In this plot, the 'p' and 'q' values can be determined as follows
Step5:  Measure the variance between the data and the values predicted by the model
You can measure whether the results of your model fit the underlying data by using the residual sum of squares (RSS) metric. A small RSS indicates that the model fits tightly to the data.
Yet another approach to validate the ARIMA model appropriateness is by performing residual analysis.
Print the results of the ARIMA model and plot the residuals. A density plot of the residual error values indicates a normal distribution centered around zero mean. Also, the residuals do not violate the assumptions of constant location and scale with most values in the range (-1,1).
Step6:  Scale predictions
Now that the model is returning the results you want to see, you can scale the model predictions back to the original scale. For this, you will remove the first order differencing and take exponent to restore the predictions back to their original scale.
The lower the root mean square error (RMSE) and the closer it is to 0, the better are the model predictions in being closer to actual values.
Step7:  Making Predictions
<ASSISTANT_TASK:>
Python Code:
test_stationary(usa_log)
Explanation: Because the p-value is larger than 0.05, the moving average is not constant over time and the null hypothesis of the Dickey-Fuller test cannot be rejected. This shows that the weekly time series is not stationary. Before you can apply ARIMA models for forecasting, you need to transform this time series into a stationary time series.
End of explanation
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(usa)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.subplot(411)
plt.plot(usa, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal, label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
Explanation: Because the p-value is now smaller than 0.05, we have necessary evidence to reject the null hypothesis of the Dickey-Fuller test. This shows that the time series is stationary at log scale. We can apply ARIMA model forcasting on this log scaled data.
Remove trend and seasonality with decomposition¶
End of explanation
usa_diff = usa - usa.shift(periods=1)
usa_diff.dropna(inplace=True)
plt.plot(usa_diff)
test_stationary(usa_diff)
Explanation: Now that you have stationarized your time series, you could go on and model residuals (fit lines between values in the plot). However, as the patterns for the trend and seasonality information extracted from the series that are plotted after decomposition are still not consistent and cannot be scaled back to the original values, you cannot use this approach to create reliable forecasts.
In case your time series does show a strong and consistent seasonal trend, here is a good article that describes the use of the SARIMAX model in Python.
Remove trend and seasonality with differencing
End of explanation
#ACF and PACF plots
plt.figure(figsize=(15, 6))
def plot_acf_pacf(ts):
    lag_acf = acf(ts, nlags=10)
    lag_pacf = pacf(ts, nlags=10, method='ols')
    #Plot ACF: 
    plt.subplot(121) 
    plt.plot(lag_acf)
    plt.axhline(y=0,linestyle='--',color='gray')
    plt.axhline(y=-7.96/np.sqrt(len(ts)),linestyle='--',color='gray')
    plt.axhline(y=7.96/np.sqrt(len(ts)),linestyle='--',color='gray')
    plt.title('Autocorrelation Function')
    #Plot PACF:
    plt.subplot(122)
    plt.plot(lag_pacf)
    plt.axhline(y=0,linestyle='--',color='gray')
    plt.axhline(y=-7.96/np.sqrt(len(ts)),linestyle='--',color='gray')
    plt.axhline(y=7.96/np.sqrt(len(ts)),linestyle='--',color='gray')
    plt.title('Partial Autocorrelation Function')
    plt.tight_layout()
    
usa_diff = usa - usa.shift(periods=1)
usa_diff.dropna(inplace=True)
plot_acf_pacf(usa_diff)
Explanation: Find optimal parameters and build an ARIMA (auto regressive integrated moving average) model
To apply an ARIMA model to your time series, you need to find optimal values for the following three model parameters (p,d,q):
 - The number of autoregressive (AR) terms (p): AR terms are just lags of the dependent variable, the euro rate, in this case. So, if p=2, it means that predictors of x(t) will be x(t-1) and x(t-2).
 - The number of moving average (MA) terms (q): MA terms are lagged forecast errors in the prediction equation. For instance, if q=2, the predictors for x(t) will be e(t-1) and e(t-2) where e(i) is the difference between the moving average at i-th instant and the actual value.
 - The number of differences (d): These are the number of non-seasonal differences. In your case, d=1, as you are modeling using the first order differenced time series.
There are two ways to determine the number of AR and MA terms. The first is to use the arma_order_select_ic function in Python. The second uses plots of the autocorrelation function (ACF) and partial autocorrelation function (PACF).
This article describes in detail the purpose of the ACF and PACF plots.
End of explanation
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(usa, order=(2, 1, 0))  
results_ARIMA = model.fit(disp=-1)  
plt.plot(usa_diff, label = "Original")
plt.plot(results_ARIMA.fittedvalues, color='red', label = "fitted")
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues - usa_diff)**2))
plt.legend()
Explanation: In this plot, the 'p' and 'q' values can be determined as follows:
 - p: The lag value where the PACF cuts off (drops to 0) for the first time. If you look closely, p=2.
 - q: The lag value where the ACF chart crosses the upper confidence interval for the first time. If you look closely, q=0.
This means that the optimal values for the ARIMA(p,d,q) model are (2,1,0).
If your assessment of the ACF and PACF plots differs from the values suggested by the arma_order_select_ic function, you should plug in different values for the p and q terms and use the model fit results to study the AIC values and proceed with the model with a lower AIC value
Run the next code cell to plot the ARIMA model using the values (2,1,0):
End of explanation
print(results_ARIMA.summary())
residuals = pd.DataFrame(results_ARIMA.resid)
residuals.plot(kind='kde')
print(residuals.describe())
Explanation: Measure the variance between the data and the values predicted by the model
You can measure whether the results of your model fit the underlying data by using the residual sum of squares (RSS) metric. A small RSS indicates that the model fits tightly to the data.
Yet another approach to validate the ARIMA model appropriateness is by performing residual analysis.
Print the results of the ARIMA model and plot the residuals. A density plot of the residual error values indicates a normal distribution centered around zero mean. Also, the residuals do not violate the assumptions of constant location and scale with most values in the range (-1,1).
End of explanation
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
predictions_ARIMA_diff.head()
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
predictions_ARIMA = pd.Series(usa.ix[0], index=usa.index)
predictions_ARIMA = predictions_ARIMA.add(predictions_ARIMA_diff_cumsum, fill_value=0)
plt.plot(usa, label = "Actual")
plt.plot(predictions_ARIMA, label = "ARIMA fitted")
plt.legend()
Explanation: Scale predictions
Now that the model is returning the results you want to see, you can scale the model predictions back to the original scale. For this, you will remove the first order differencing and take exponent to restore the predictions back to their original scale.
The lower the root mean square error (RMSE) and the closer it is to 0, the better are the model predictions in being closer to actual values.
End of explanation
from sklearn.metrics import mean_squared_error
size = int(len(usa) - 5)
train, test = usa[0:size], usa[size:len(usa)]
history = [x for x in train]
predictions = list()
print('Printing Predicted vs Expected Values...')
print('\n')
for t in range(len(test)):
    model = ARIMA(history, order=(2,1,0))
    model_fit = model.fit(disp=0)
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(float(yhat))
    obs = test[t]
    history.append(obs)
    print('predicted=%f, expected=%f' % (yhat, obs))
error = mean_squared_error(test, predictions)
print('\n')
print('Printing Mean Squared Error of Predictions...')
print('Test MSE: %.6f' % error)
predictions_series = pd.Series(predictions, index = test.index)
#rolling one-step out-of-sample
fig, ax = plt.subplots()
ax.set(title='GDP Per Capita forecasting', xlabel='Date', ylabel='Per Capita GDP')
ax.plot(usa, marker = '.', color = "red", label='observed')
ax.plot(predictions_series, color = "green", label='forecast')
legend = ax.legend(loc='upper left')
legend.get_frame().set_facecolor('w')
Explanation: Making Predictions
End of explanation
<END_TASK>