# ARMA Models

## Intro to time series and stationarity


#### Exploration


```python
# Import modules
import matplotlib.pyplot as plt
import pandas as pd

# Load in the time series
candy = pd.read_csv('candy_production.csv', 
            index_col='date',
            parse_dates=True)

# Plot and show the time series on axis ax1
fig, ax1 = plt.subplots()
candy.plot(ax=ax1)
plt.show()
```

#### Train-test split


```python
# Split the data into a train and test set
candy_train = candy.loc[:'2006']
candy_test = candy.loc['2007':]

# Create an axis
fig, ax = plt.subplots()

# Plot the train and test sets on the axis ax
candy_train.plot(ax=ax)
candy_test.plot(ax=ax)
plt.show()
```

### Making time series stationary


#### Testing for stationarity wit Augmented Dicky-Fuller



```python
# Import augmented dicky-fuller test function
from statsmodels.tsa.stattools import adfuller

# Run test
result = adfuller(earthquake['earthquakes_per_year'])

# Print test statistic
print(result[0])

# Print p-value
print(result[1])

# Print critical values
print(result[4]) 
```

#### Taking the difference
* .diff() applyed twice

#### Other tranforms

Differencing should be the first transform you try to make a time series stationary. But sometimes it isn't the best option.

A classic way of transforming stock time series is the log-return of the series. This is calculated as follows:

\begin{equation}
log\_return(y_t) = log\frac{y_{t}}{y_{t-1}}
\end{equation}

The Amazon stock time series has already been loaded for you as amazon. You can calculate the log-return of this DataFrame by substituting:

 * y_t -> amazon
 * y_(t-1) -> amazon.shift(1)
 * log -> np.log()

In this exercise you will compare the log-return transform and the first order difference of the Amazon stock time series to find which is better for making the time series stationary.


```python
# Calculate the first difference and drop the nans
amazon_diff = amazon.diff()
amazon_diff = amazon_diff.dropna()

# Run test and print
result_diff = adfuller(amazon_diff['close'])
print(result_diff)

# Calculate log-return and drop nans
amazon_log = np.log(amazon.div(amazon.shift(1)))
amazon_log = amazon_log.dropna()

# Run test and print
result_log = adfuller(amazon_log['close'])
print(result_log)
```

## Structure of ARIMA models

#### AR Models

#### MA Models

#### ARMA Models

#### ARMAX Models

* X - Exogenous variable
* ARMAX = ARMA + linear regression

#### Generating ARMA data



```python
# Import data generation function and set random seed
from statsmodels.tsa.arima_process import arma_generate_sample
np.random.seed(3)

# Set coefficients
ar_coefs = [1, 0.2] # Negative coeficients
ma_coefs = [1, 0.3, 0.4]

# Generate data
y = arma_generate_sample(ar_coefs, ma_coefs, nsample=100, scale=0.5)

plt.plot(y)
plt.ylabel(r'$y_t$')
plt.xlabel(r'$t$')
plt.show()
```

#### Fitting Prelude



```python
# Import the ARMA model
from statsmodels.tsa.arima_model import ARMA

# Instantiate the ARMA model 
model = ARMA(y, order=(1,1))

# Instantiate the ARMAX model
model = ARMA(S1, order = (2,1), exog = S2)

# Fit the model
results = model.fit()

# Print summary
print(results.summary())
```

# Fitting the Future

#### Generating one-step-ahead predictions



```python
# Generate predictions
one_step_forecast = results.get_prediction(start=-30)

# Extract prediction mean
mean_forecast = one_step_forecast.predicted_mean

# Get confidence intervals of predictions
confidence_intervals = one_step_forecast.conf_int()

# Select lower and upper confidence limits
lower_limits = confidence_intervals.loc[:,'lower close']
upper_limits = confidence_intervals.loc[:,'upper close']

# Print best estimate predictions
print(mean_forecast)
```

#### Plotting one-step-ahead predictions



```python
# plot the amazon data
plt.plot(amazon.index, amazon, label='observed')

# plot your mean predictions
plt.plot(mean_forecast.index, mean_forecast, color='r', label='forecast')

# shade the area between your confidence limits
plt.fill_between(lower_limits.index, lower_limits, 
               upper_limits, color='pink')

# set labels, legends and show plot
plt.xlabel('Date')
plt.ylabel('Amazon Stock Price - Close USD')
plt.legend()
plt.show()
```

#### Generating dynamic forecasts



```python
# Generate predictions
dynamic_forecast = results.get_prediction(start = -30, dynamic=True)

# Extract prediction mean
mean_forecast = dynamic_forecast.predicted_mean

# Get confidence intervals of predictions
confidence_intervals = dynamic_forecast.conf_int()

# Select lower and upper confidence limits
lower_limits = confidence_intervals.loc[:,'lower close']
upper_limits = confidence_intervals.loc[:,'upper close']

# Print best estimate predictions
print(mean_forecast)
```

#### Differencing and fitting ARMA


```python
# Take the first difference of the data
amazon_diff = amazon.diff().dropna()

# Create ARMA(2,2) model
arma = SARIMAX(amazon_diff, order = (2,0,2))

# Fit model
arma_results = arma.fit()

# Print fit summary
print(arma_results.summary())
```

#### Unrolling ARMA forecast


```python
# Make arma forecast of next 10 differences
arma_diff_forecast = arma_results.get_forecast(steps=10).predicted_mean

# Integrate the difference forecast
arma_int_forecast = np.cumsum(arma_diff_forecast)

# Make absolute value forecast
arma_value_forecast = arma_int_forecast + amazon.iloc[-1,0]

# Print forecast
print(arma_value_forecast)
```

#### Using ARIMA


```python
# Create ARIMA(2,1,2) model
arima = SARIMAX(amazon, order = (2,1,2))

# Fit ARIMA model
arima_results = arima.fit()

# Make ARIMA forecast of next 10 values
arima_value_forecast = arima_results.get_forecast(steps=10).predicted_mean

# Print forecast
print(arima_value_forecast)
```

# The Best of the Best Models

## ACF & PCF (auto correlation and parcial auto correlation)

* IF (ACF tails off) and (PACF cuts off after lag p): AR(p)
* IF (ACF cuts off after q) and (APCF tails off): MA(q)
* IF both tail: ARMA
* IF (ACF values are high) and (tail slowly): non stationary
* IF (ACF lag-1 = too negative): you took the difference too many times

## AIC & BIC

* AIC: Akaike information criteria
* BIC: Bayesian information criteria

Lower AIC and BIC indicate better models. They penalize excess of parameters.

BIC focusses more on restricting model size and is better for choosing a explanatory model.

AIC is better at choosing predictive models.

#### Searching over model order



```python
# Create empty list to store search results
order_aic_bic=[]

# Loop over p values from 0-2
for p in range(3):
  # Loop over q values from 0-2
    for q in range(3):
      	# create and fit ARMA(p,q) model
        model = SARIMAX(df, order=(p,0,q))
        results = model.fit()
        
        # Append order and results tuple
        order_aic_bic.append((p,q,results.aic, results.bic)),
        
# Construct DataFrame from order_aic_bic
order_df = pd.DataFrame(order_aic_bic, 
                        columns=['p', 'q', 'AIC', 'BIC'])

# Print order_df in order of increasing AIC
print(order_df.sort_values('AIC'))

# Print order_df in order of increasing BIC
print(order_df.sort_values('BIC'))
```

## Model Diagnostics

#### Mean absolute error



```python
# Fit model
model = SARIMAX(earthquake, order=(1,0,1))
results = model.fit()

# Calculate the mean absolute error from residuals
mae = np.mean(np.abs(results.resid))

# Print mean absolute error
print(mae)

# Make plot of time series for comparison
earthquake.plot()
plt.show()
```

#### Diagnostic summary statistics

It is important to know when you need to go back to the drawing board in model design. In this exercise you will use the residual test statistics in the results summary to decide whether a model is a good fit to a time series.

Here is a reminder of the tests in the model summary:

| Test | Null hypothesis  | P-value name  | IF|
|---|---|---|---|
|  Ljung-Box |  There are no correlations in the residual |  Prob(Q) | p<0.05 : correlated|
|  Jarque-Bera |  The residuals are normally distributed | Prob(JB)  | p<0.05 : non normal|

We reject the null hypothesis if the p-value < 0.05.
		


```python
# Create and fit model
model1 = SARIMAX(df, order=(3,0,1))
results1 = model1.fit()

# Print summary
print(results1.summary())
```

#### Plot diagnostics

|Test	| Good fit |
|---|---|
|Standardized residual |	There are no obvious patterns in the residuals|
|Histogram plus kde estimate | The KDE curve should be very similar to the normal distribution|
|Normal Q-Q	 | Most of the data points should lie on the straight line|
|Correlogram	| 95% of correlations for lag greater than zero should not be significant|


```python
# Create and fit model
model = SARIMAX(df, order=(1,1,1))
results = model.fit()

# Create the 4 diagostics plots
results.plot_diagnostics()
plt.show()
```

### Box-Jenkins method

From Raw data -> production model

* Identification
* Estimation
* Model siagnostics

#### Identification

* Is the time series stationary?
* What differencing will make it stationary?
* What transforms will make it stationary?
* What values of P and q are most promissing?
    * Plot the time series
    * Use audmented Dicky-fuller test
    * Use transforms .diff(), np.log(), np.sqrt()
    * Plot ACF/PACF and pick a model order.

#### Estimation

* Use data to train model coefficients
* Choose between models using AIC and BIC

#### Model Diagnostics

* Are residuals uncorrelated?
* Are residuals normally distributed?
    * results.plot_diagnostics()
    * results.summary()

#### Is model ok?
* Back to identification
* Send to production

#### Identification


```python
# Plot time series
savings.plot()
plt.show()

# Run Dicky-Fuller test
result = adfuller(savings.savings)

# Print test statistic
print(result[0])

# Print p-value
print(result[1])

# Create figure
fig, (ax1, ax2) = plt.subplots(2,1, figsize=(12,8))
 
# Plot the ACF of savings on ax1
plot_acf(savings, lags = 10, zero = False, ax = ax1)

# Plot the PACF of savings on ax2
ax2 = plot_pacf(savings, lags = 10,  zero = False, ax = ax2)

plt.show()
```

#### Estimation


```python
# Loop over p values from 0-3
for p in range(4):
  
  # Loop over q values from 0-3
    for q in range(4):
      try:
        # Create and fit ARMA(p,q) model
        model = SARIMAX(savings, order=(p,0,q), trend='c')
        results = model.fit()
        
        # Print p, q, AIC, BIC
        print(p, q, results.aic, results.bic)
        
      except:
        print(p, q, None, None)
```

#### Diagnostics


```python
# Create and fit model
model = SARIMAX(savings, order = (1, 0, 2), trend = 'c')
results = model.fit()

# Create the 4 diagostics plots
results.plot_diagnostics()
plt.show()

# Print summary
print(results.summary())
```

# Seasonal ARIMA Models (SARIMA)

#### Seasonal decompose


```python

# Import seasonal decompose
from statsmodels.tsa.seasonal import seasonal_decompose

# Perform additive decomposition
decomp = seasonal_decompose(milk_production['pounds_per_cow'], 
                            period=12)

# Plot decomposition
decomp.plot()
plt.show()
```

#### Seasonal ACF and PACF



```python
# Create figure and subplot
fig, ax1 = plt.subplots()

# Plot the ACF on ax1
plot_acf(water['water_consumers'], lags = 25, zero=False,  ax=ax1)

# Show figure
plt.show()

# Subtract the rolling mean
water_2 = water - water.rolling(15).mean()

# Drop the NaN values
water_2 = water_2.dropna()

# Create figure and subplots
fig, ax1 = plt.subplots()

# Plot the ACF
plot_acf(water_2['water_consumers'], lags=25, zero=False, ax=ax1)

# Show figure
plt.show()

# Components higher than the confidence interval means seasonal components.
```

## SARIMA models

To fit a seasonal ARIMA, we fit 2 models:
* One for the seasonal part.
* One for the non seasonal part.

\begin{equation}
SARIMA(p,d,q)(P,D,Q)_{S}
\end{equation}
S: number of time steps per cycle

### Seasonal differencing
To make a time searies stationary, subtract a value from a season ago from the time series.

To find the seasonal order, plot the ACF and PACF of the differeced time series at multiple season steps.

#### Fitting SARIMA models


```python
# Create a SARIMAX model
model = SARIMAX(df1, order=(1,0,0), seasonal_order=(1,1,0,7))

# Fit the model
results = model.fit()

# Print the results summary
print(results.summary())

# Create a SARIMAX model
model = SARIMAX(df2, order = (2,1,1), seasonal_order = (1,0,0,4))

# Fit the model
results = model.fit()

# Print the results summary
print(results.summary())

# Create a SARIMAX model
model = SARIMAX(df3,  order = (1,1,0), seasonal_order = (0,1,1,12))

# Fit the model
results = model.fit()

# Print the results summary
print(results.summary())
```

#### Choosing SARIMA order



```python
# Take the first and seasonal differences and drop NaNs
aus_employment_diff = aus_employment.diff(1).diff(12).dropna()

# Create the figure 
fig, (ax1, ax2) = plt.subplots(2,1,figsize=(8,6))

# Plot the ACF on ax1
plot_acf(aus_employment_diff, lags = 11, zero = False, ax = ax1)

# Plot the PACF on ax2
plot_pacf(aus_employment_diff, lags = 11, zero = False, ax = ax2)

plt.show()

# Make list of lags
lags = [12, 24, 36, 48, 60]

# Create the figure 
fig, (ax1, ax2) = plt.subplots(2,1,figsize=(8,6))

# Plot the ACF on ax1
plot_acf(aus_employment_diff, lags = lags, zero = False, ax = ax1)

# Plot the PACF on ax2
plot_pacf(aus_employment_diff, lags = lags, zero = False, ax = ax2)

plt.show()

```

#### SARIMA vs ARIMA forecasts

# Create ARIMA mean forecast
arima_pred = arima_results.get_prediction(start = -25)
arima_mean = arima_pred.predicted_mean

# Create SARIMA mean forecast
sarima_pred = sarima_results.get_prediction(start = -25)
sarima_mean = sarima_pred.predicted_mean

# Plot mean ARIMA and SARIMA predictions and observed
plt.plot(dates, sarima_mean, label='SARIMA')
plt.plot(dates, arima_mean, label='ARIMA')
plt.plot(wisconsin_test, label='observed')
plt.legend()
plt.show()

# Automation and saving to production

#### Automated model selection



```python
# Import pmdarima as pm
import pmdarima as pm

# Create auto_arima model
model1 = pm.auto_arima(df1,
                      seasonal=True, m=7,
                      d=0, D=1, 
                 	  max_p=2, max_q=2,
                      trace=True,
                      error_action='ignore',
                      suppress_warnings=True)
                       
# Print model summary
print(model1.summary())
```


```python
# Create model
model2 = pm.auto_arima(df2,
                      seasonal=False, 
                      d=1, 
                      trend='c',
                 	  max_p=2, max_q=2,
                      trace=True,
                      error_action='ignore',
                      suppress_warnings=True) 

# Print model summary
print(model2.summary())
```


```python
# Create model for SARIMAX(p,1,q)(P,1,Q)7
model3 = pm.auto_arima(df3,
                      seasonal=True, m=7,
                      d=1, D=1, 
                      start_p=1, start_q=1,
                      max_p=1, max_q=1,
                      max_P=1, max_Q=1,
                      trace=True,
                      error_action='ignore',
                      suppress_warnings=True) 

# Print model summary
print(model3.summary())
```

#### Saving and updating models



```python
# Import joblib
import joblib

# Set model name
filename = "candy_model.pkl"

# Pickle it
joblib.dump(model,filename)

# Load the model back in
loaded_model = joblib.load(filename)

# Update the model
loaded_model.update(df_new)
```

# Box-Jenkins with seasonal data

D should be 0 or 1

d + D should be 0-2 

Take the log of the data if the relations between seasonal and non seasonal componentes are multiplicative.



# Final Example

#### Sarima model diagnostics


```python
# Import model class
from statsmodels.tsa.statespace.sarimax import SARIMAX

# Create model object
model = SARIMAX(co2, 
                order=(1,1,1), 
                seasonal_order=(0,1,1,12), 
                trend='c')
# Fit model
results = model.fit()

# Plot common diagnostics
results.plot_diagnostics()
plt.show()
```


```python
# Create forecast object
forecast_object = results.get_forecast(steps = 136)

# Extract predicted mean attribute
mean = forecast_object.predicted_mean

# Calculate the confidence intervals
conf_int = forecast_object.conf_int()

# Extract the forecast dates
dates = mean.index

plt.figure()

# Plot past CO2 levels
plt.plot(co2.index, co2, label='past')

# Plot the prediction means as line
plt.plot(dates, mean, label='predicted')

# Shade between the confidence intervals
plt.fill_between(dates, conf_int['lower CO2_ppm'], conf_int['upper CO2_ppm'], alpha=0.2)

# Plot legend and show figure
plt.legend()
plt.show()

# Print last predicted mean
print(mean.iloc[-1])

# Print last confidence interval
print(conf_int.iloc[-1])
```

b9178060035cff63949c836be03bee97dcfb8297
ipynb
Jupyter Notebook
ARIMA Models in Python.ipynb
enerdb/Data-Camp-Courses
871501858afb696a12627c310870cb0d25bf9acc
MIT
2022-03-12T23:33:07.000Z
2022-03-12T23:33:07.000Z
ARIMA Models in Python.ipynb
enerdb/Data-Camp-Courses
871501858afb696a12627c310870cb0d25bf9acc
MIT
ARIMA Models in Python.ipynb
enerdb/Data-Camp-Courses
871501858afb696a12627c310870cb0d25bf9acc
MIT
Qwen/Qwen-72B
1. YES
2. YES