Time Series Analysis and SARIMA Model for Stock Price Prediction

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import yfinance as yf
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller, acf, pacf
from pmdarima import auto_arima
import seaborn as sns


# Download Apple stock data
ticker = "AAPL"
start_date = "2022-01-01"
end_date = "2024-01-23"
data = yf.download(ticker, start=start_date, end=end_date)
df = data.copy()

# Use closing prices for our analysis
ts = data['Close'].dropna()

# Plot the time series
plt.figure(figsize=(12,6))
plt.plot(ts)
plt.title('Apple Stock Closing Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.show()

# Perform seasonal decomposition
decomposition = seasonal_decompose(ts, model='additive', period=12)  # 252 trading days in a year
plt.rcParams.update({'figure.dpi':200})

# Plot the decomposition
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(12, 16))
decomposition.observed.plot(ax=ax1)
ax1.set_title('Observed')
decomposition.trend.plot(ax=ax2, c='r')
ax2.set_title('Trend')
decomposition.seasonal.plot(ax=ax3, c='k')
ax3.set_title('Seasonal')
decomposition.resid.plot(ax=ax4, c='g')
ax4.set_title('Residual')
plt.tight_layout()
plt.show()

from pmdarima.arima.utils import nsdiffs
D = nsdiffs(ts,
            m=12,  # commonly requires knowledge of dataset
            max_D=12,
            )  # -> 0
D

# Create Training and Test
train = ts[:int(len(df)*0.8)]
test = ts[int(len(df)*0.8):]

def grid_search_sarima(train, test, p_range, d_range, q_range, P_range, D_range, Q_range, m_range):
    best_aic = float('inf')
    best_mape = float('inf')
    best_order = None
    for p in p_range:
        for d in d_range:
            for q in q_range:
                for P in P_range:
                    for D in D_range:
                        for Q in Q_range:
                            for m in m_range:
                                try:
                                    model = SARIMAX(train.values, order=(p,d,q), seasonal_order=(P,D,Q,m))
                                    results = model.fit()
                                    fc_series = pd.Series(results.forecast(steps=len(test)), index=test.index)  # 95% conf
                                    test_metrics = forecast_accuracy(fc_series.values, test.values)
                                    # if results.aic < best_aic:
                                    #     best_aic = results.aic
                                    #     best_order = (p,d,q)
                                    # print(f"(p,d,q)=({p},{d},{q}) | (P,D,Q,m)=({P},{D},{Q},{m})", test_metrics['mape'])
                                    if test_metrics['mape'] < best_mape:
                                        best_mape = test_metrics['mape']
                                        best_order = (p,d,q)
                                        best_sorder = (P,D,Q,m)
                                        print("temp best:" + f"(p,d,q)=({p},{d},{q}) | (P,D,Q,m)=({P},{D},{Q},{m})", 
                                            round(test_metrics['mape'], 4))
                                except Exception as e:
                                    print(e)
                                    continue
    return best_order, best_sorder


# Accuracy metrics
def forecast_accuracy(forecast, actual):
    mape = np.mean(np.abs(forecast - actual)/np.abs(actual))    # MAPE
    me = np.mean(forecast - actual)                             # ME
    mae = np.mean(np.abs(forecast - actual))                    # MAE
    mpe = np.mean((forecast - actual)/actual)                   # MPE
    rmse = np.mean((forecast - actual)**2)**.5                  # RMSE
    corr = np.corrcoef(forecast, actual)[0,1]                   # corr
    mins = np.amin(np.hstack([forecast[:,None], 
                              actual[:,None]]), axis=1)
    maxs = np.amax(np.hstack([forecast[:,None], 
                              actual[:,None]]), axis=1)
    minmax = 1 - np.mean(mins/maxs)                             # minmax
    acf1 = acf(forecast-test)[1]                                # ACF1
    return({'mape':mape, 'me':me, 'mae': mae, 
            'mpe': mpe, 'rmse':rmse, 'acf1':acf1, 
            'corr':corr, 'minmax':minmax})

# Determine optimal SARIMA parameters
# model = auto_arima(ts, seasonal=True, m=12,
#                     start_p=0, start_q=0, start_P=1, start_Q=1,
#                     max_p=4, max_q=4, max_P=2, max_Q=2, d=1, D=1,
#                     trace=True, error_action='ignore', suppress_warnings=True, stepwise=True, out_of_sample=len(test))

# print(model.summary())

# Custom Grid Search
# d = 1 # from 
# best_order, best_sorder = grid_search_sarima(train, test, range(0,3), [d], range(0,3), range(3), [D], range(3), [4, 8, 12])
# print(f"Best ARIMA order based on grid search: {best_order}")

# Fit the SARIMA model
# sarima_model = SARIMAX(ts, order=model.order, seasonal_order=model.seasonal_order) (2,1,2) | (P,D,Q,m)=(2,1,1,
# sarima_model = SARIMAX(train, order=best_order, seasonal_order=best_sorder) #(2,1,2) | (P,D,Q,m)=(2,1,1,8)
# print(best_order, best_sorder)
# sarima_model = SARIMAX(train, order=(2,1,2), seasonal_order=(6,0,6,8))
sarima_model = SARIMAX(train, order=(2,1,2), seasonal_order=(2,2,1,8))
results = sarima_model.fit()

print(results.summary())

# sarima_model = SARIMAX(train, order=(2,1,2), seasonal_order=(2,2,1,8))
# results = sarima_model.fit()


# Forecast
forecast_steps = len(test) # Forecast for one year
forecast = results.get_forecast(steps=forecast_steps)
forecast_ci = forecast.conf_int(alpha=0.25)
# forecast_index = pd.date_range(start=test.index[0]+pd.Timedelta(24, 'hours'), periods=len(forecast_ci))
forecast_index = test.index

# Plot the forecast
plt.figure(figsize=(12, 6))
plt.plot(train[-200:], label='training')
# plt.plot(ts.index, ts, label='Observed')
plt.plot(test, label='actual')
plt.plot(forecast_index, forecast.predicted_mean + decomposition.seasonal[test.index].values*10, color='r', label='forecast')
plt.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='pink', alpha=0.4)
plt.title(f'SARIMA Forecast of {ticker} Stock Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.legend()
plt.show()

# Model diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.show()