data.head() # that's better

data.shape

data.tail() # need to deal with the last row

data.tail() # that's better

data.dtypes # datatypes shouldn't be "object"

data.dtypes # no longer "object"

data.head() # It looks like the Riders header is over the Month header.  My guess is the Month column is the index for the dataframe.

data.describe()

# Plot the riders over time
data.plot(figsize=(12,6))

# Split the results into its Systematic (Seasonal and Non-Seasonal) and Non-Systematic components by using the seasonal_decompose function
results = seasonal_decompose(data['Riders']);
results.plot();    

# EXTRA: Additive model
resultsAdd = seasonal_decompose(data['Riders'], model = 'additive');
resultsAdd.plot(); # This looks the same as above, so apparently that is the default for the seasonal_decompose function  

# EXTRA: Multiplicitive model
resultsMult = seasonal_decompose(data['Riders'], model = 'multiplicative');
resultsMult.plot(); # Actually, this looks a lot more similar to the Additive model than I was expecting.

# Do an ADF Test and output results
adfuller_test(data['Riders'])

# Plot the data, with indicators for the years
ax = data.plot(figsize=(12,6))

xcoords = ['1960-01-01', '1961-01-01','1962-01-01', '1963-01-01', '1964-01-01', '1965-01-01',
          '1966-01-01', '1967-01-01', '1968-01-01', '1969-01-01']
for xc in xcoords:
    plt.axvline(x=xc, color='black', linestyle='--')

# Plotting the Autocorrelation and Partial Autocorrelation.
# This page has some info about these two terms: https://statisticsbyjim.com/time-series/autocorrelation-partial-autocorrelation/
# Autocorrelation is the correlation between two values in a time series.
# The partial autocorrelation function is similar to the ACF, 
#    except that it displays only the correlation between two observations that the shorter lags between those observations do not explain.
testing = data['Riders'].dropna()
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(testing, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(testing, lags=40, ax=ax2)
# Since the points in the Autocorrelation plot are tapering off slowly, this indicates that there is a trend to the data.
# However, the points are going pretty steadily down instead of waving, which indicates not much seasonality.
# I thought the plots above looked pretty seasonal though.
# Since the first 2 points in the partial autocorrelation plot are above the shaded part, 
#    this indicates fitting a 2nd order autogressive model

# This page explains about autogressive models: https://sjmiller8182.github.io/LearningTS/build/autoregressive_models.html
# The interpretation of AR(2) is 
#    "the value at time t is a linear combination of two values at the two previous time periods plus a random noise component".

# Check to make sure the shift was added correctly
data.head(20)

# Check if the data is stationary when looking at the change between years
adfuller_test(data['Seasonal_1st_diff'].dropna())
# It still says that the data is not stationary, but the p value is closer to .05

# Plot the season difference
data['Seasonal_1st_diff'].plot()

# Add the year lines again
ax = data['Seasonal_1st_diff'].plot(figsize=(12,6), linewidth=3, fontsize=20)

xcoords = ['1960-01-01', '1961-01-01','1962-01-01', '1963-01-01', '1964-01-01', '1965-01-01',
          '1966-01-01', '1967-01-01', '1968-01-01', '1969-01-01']
for xc in xcoords:
    plt.axvline(x=xc, color='black', linestyle='--')

# Plot the autocorrelation
# This page explains what Panda's autocorrelation_plot function does: https://blog.finxter.com/pandas-plotting-autocorrelation/
autocorrelation_plot(data['Riders'])
plt.show()
# Since it has some values above the dotted lines at +/- 0.25, it indicates there is a trend.

# EXTRA: What does the autocorrelation_plot look like for the seasonal diff?
testing2 = data['Seasonal_1st_diff'].dropna()
autocorrelation_plot(testing2)
plt.show()
# Still a drop at the beginning, but a little bit more random looking

# EXTRA: What does the original autocorrelation function look like?
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(testing2, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(testing2, lags=40, ax=ax2)
# The trend is still there

# Try drawing a lag plot
# Here is an explantion of the lag_plot in Python: https://pythontic.com/pandas/plotting/lag%20plot
pd.plotting.lag_plot(data['Riders'])

# EXTRA: Try a lag plot with the seasonal shift data
pd.plotting.lag_plot(data['Seasonal_1st_diff'])
# It looks like a less tight correlation than the original data has

# Determine what p, d, and q should be for the arima model
stepwise_fit = auto_arima(data['Riders'], trace=True)
# It looks like p = 1, d = 1, q = 3

# One of the articles that I had looked at before for general information about the SARIMAX algorithm talked about p, d, and q: 
#    https://medium.com/swlh/a-brief-introduction-to-arima-and-sarima-modeling-in-python-87a58d375def
# p = auto regressive (AR) order: determines how many previous data points will be used
# d: how many times the data is to be differenced
# q = moving average (MA) order: determines how many previous data points will be used

# Look at the metrics for the chosen ARIMA model
stepwise_fit.summary()

# Make the SARIMAX model
# NOTE: It makes sense to use a SARIMAX model here because we noticed seasonality, but even with the seasonality
model = sm.tsa.statespace.SARIMAX(data['Riders'], order=(1,1,3), seasonal_order = (1,1,3,12))
results = model.fit()

# Plot predicted vs actual for the last couple years
data['forecast'] = results.predict(start=90, end=114, dynamic=True)
data[['Riders', 'forecast']].plot(figsize=(12,8), linewidth=3, fontsize=20)
# The predicted values seem to be much higher than the actual values.  
# It does appear to be a drop of after 1967, so that makes some sense.

# See what the last dates are
future_dates_df.tail()
# Looks like the forecast dates go to May 1971

# Combine the datasets
future_data = pd.concat([data, future_dates_df])
future_data

# Plotting the original data with the forecasted data
future_data['forecast'] =  results.predict(start=113, end=128, dynamic=True)
future_data[['Riders', 'forecast']].plot(figsize=(12,8), linewidth=3, fontsize=20)

# What is the error for the original parameter set
rse

print(f'Final parameters: p1 = {best_p1}, q1 = {best_q1}, p2 = {best_p2}, q2 = {best_q2}. Best RSE = {best_rse}')

# Plotting the two forecasts
model2 = sm.tsa.statespace.SARIMAX(data['Riders'], order=(4,1,1), seasonal_order = (0,1,0,12))
results2 = model2.fit()
data['forecast2'] = results2.predict(start=90, end=114, dynamic=True)
data[['Riders', 'forecast', 'forecast2']].plot(figsize=(12,8), linewidth=3, fontsize=20)
# Forecast2 is a little bit less extreme than Forecast

# Plot the future forcasts from forecast and forecast2
future_data['forecast2'] =  results2.predict(start=113, end=128, dynamic=True)
future_data[['Riders', 'forecast', 'forecast2']].plot(figsize=(12,8), linewidth=3, fontsize=20)
# I don't like how the forecast2 prediction seems to be decreasing, even though there has been a clear upward trend from 1960 to 1967.