<jupyter_start><jupyter_text>Autoregressive model using OLSIn this notebook we will use Ordinary Least Square (OLS) to fit a linear model to time series data. OLS is a supervised learning approach and is not designed specifically with time series in mind.  However, it is easy to adapt and use with autoregressive data.  There are two caveats     1. Unlike purpose built time series models such as ARIMA we need to preprocess the data beforehand.2. Forecasting **h-steps** ahead is more involved than other methods. However, OLS itself is very simple to use and the procedures we develop to solve 1 and 2 are exactly what we need to work with Feedforward Neural Networks and Recurrent Neural Networks which we will cover later.---**LEARNING OBJECTIVES*** Mechanics of fitting time series data to an OLS autoregressive model* Generate h-step forecasts using an iterative approach* Generate h-step forecast using a direct modelling approach---  1. Python dependenciesIt is recommended that you use the forecasting course conda environment provided for this work.  The implementation of Ordinary Least Squares Regression that we are going to use is in `statsmodels`.  You should be using at least version `0.11.1`<jupyter_code>import statsmodels.api as sm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sm.__version__<jupyter_output><empty_output><jupyter_text>The `statsmodels` package has a class for OLS that we are going to use.```python statsmodels.regression.linear_model.OLS ```<jupyter_code>from statsmodels.regression.linear_model import OLS<jupyter_output><empty_output><jupyter_text>2. AR OLS: The forecasting process1. Select $l$ the number of autoregressive lags and forecast horizon $h$2. Preprocess the data into tabular (supervised learning) form [[$lag_1, lag_2, ... lag_l$], [$y_t$]]3. Fit the OLS model to the tabular data4. Iteratively forecast 1-step ahead gradually replacing ground truth observations with predictions. 2.1 Synthetic data without noiseGiven the extra complexities of forecasting using OLS we will use simple synthetic data before exploring real healthcare data. The synthetic data we wil use is a cosine.<jupyter_code>t = np.arange(200)
ts_data = np.cos(0.2 * t)
plt.plot(ts_data);
ts_data<jupyter_output><empty_output><jupyter_text>2.2. Preprocess the time series into tabular autoregressive formAn autoregressive model consists of $l$ lags of the time series. An easy way to think about the form of the data for a autoregressive OLS model is as a table of variables.  The first $l$ columns are the lags (the independent predictor variables) and the final column is $y$ at time $t$ ($y_t$) that is, the target/dependent variable.  We there need to manipulate the time series so that is now in that format.  More precisely for each row we need: **A vector presenting the lags at time t*** $X_t = $ [$lag_{t-l}, ... lag_{t-2}, lag_{t-1}$]**A scalar value representing y at time t:*** $y_t$For training we need a vector of rows ($X_t$) and vector of target $y_t$. e.g.```pythonX_train = [X_1, X_2, X_3, ...,  X_t]y_train = [y_1, y_2, y_3, ..., y_t]```---The function `sliding_window` illustrates how to preprocess time series data into tabular form  in python.<jupyter_code>def sliding_window(train, window_size=2, horizon=1):
    '''
    sliding window.
    
    Parameters:
    --------
    train: array-like
        training data for time series method
    
    window_size: int, optional (default=2)
        lookback - how much lagged data to include.
        
    horizon: int, optional (default=1)
        number of observations ahead to predict
            
    Returns:
        array-like, array-like
    
        preprocessed X, preprocessed Y
    '''
    tabular_X = []
    tabular_y = []
    
    for i in range(0, len(train) - window_size - horizon):
        X_train = train[i:window_size+i]
        y_train = train[i+window_size+horizon-1]
        tabular_X.append(X_train)
        tabular_y.append(y_train)
       
    return np.asarray(tabular_X), np.asarray(tabular_y).reshape(-1, 1)
#Let's do a train test split first
train = ts_data[:175]
test = ts_data[175:]

#preprocess time series training and test sets
X_train, y_train = sliding_window(train, window_size=2)
X_test, y_test = sliding_window(test, window_size=2)
print(X_train.shape)
print(y_train.shape)<jupyter_output>(172, 2)
(172, 1)<jupyter_text>The preprocessing included 2 lags.  So let's have a look at how the first elements in `X_train` and `y_train` compare to the first three elements of `ts_data`<jupyter_code>print(f'ts_data[:3]: {ts_data[:3]}')
print(f'X_train[0]: {X_train[0]}')
print(f'y_train[0]: {y_train[0]}')<jupyter_output>ts_data[:3]: [1.         0.98006658 0.92106099]
X_train[0]: [1.         0.98006658]
y_train[0]: [0.92106099]<jupyter_text>Just to prove this is work let's manually slide out window of lag size 2 along `ts_data` to see how it compares with the next element in `X_train` and `y_train`<jupyter_code>print(f'ts_data[1:4]: {ts_data[1:4]}')
print(f'X_train[1]: {X_train[1]}')
print(f'y_train[1]: {y_train[1]}')<jupyter_output>ts_data[1:4]: [0.98006658 0.92106099 0.82533561]
X_train[1]: [0.98006658 0.92106099]
y_train[1]: [0.82533561]<jupyter_text>Just to be absolutely clear what we are doing let's convert the numpy arrays to a **pandas.DataFrame**.  In this format you should be able to see why the an OLS modelcan be used for autoregression in a time series.<jupyter_code>tabular_form = pd.concat([pd.DataFrame(X_train), pd.DataFrame(y_train)], axis=1)
columns = [f'lag_{i}' for i in range(len(X_train[0]), 0, -1)]
columns.append('y_t')
tabular_form.columns = columns
tabular_form.head()<jupyter_output><empty_output><jupyter_text>2.3 Fit the preprocessed data to the OLS modelAfter preprocessing the data, fitting the data is relatively straightforward.  We create an instance of `OLS` passing in the training data and call the `.fit()` method.  The method fit returns a `RegressionResults` object that we use for prediction.  We can called the `.summary()` method to see the fitted model.<jupyter_code>X_train = sm.add_constant(X_train)
model = OLS(endog=y_train, exog=X_train).fit()

#see the regression results (not the adj. R squared!)
model.summary()<jupyter_output><empty_output><jupyter_text>2.4 Forecasting 1 step aheadTo forecast 1-step ahead we use the `RegressionResults` method `.predict(exog)`.  For example if we wanted to forecast the first y observation in the test set.<jupyter_code>X_test = sm.add_constant(X_test)
pred = model.predict(exog=X_test[0])[0]
print(f'1-step forecast: {pred}')
print(f'ground trust value: {y_test[0][0]}')<jupyter_output>1-step forecast: -0.6656134553337546
ground trust value: -0.6656134553337595<jupyter_text>An alternative way to call predict is to use `.get_prediction(exog)` to get a `PredictionResults` object.  This has the nice method `.summary_frame(alpha=0.05)` that can be used to get a dataframe with<jupyter_code>results = model.get_prediction(exog=X_test[0])
results.summary_frame(alpha=0.2)<jupyter_output><empty_output><jupyter_text>2.4 Forecast h periods ahead using the iterative method.**We have trained our `OLS` model to predict 1-step** ahead by regressing the previous two observations in the time series.  When forecasting 2 or more steps ahead we still only have two ground truth observations ($lag_1$ and $lag_2$). This means that when forecasting h-steps ahead we need to do this in a loop where we iteratively replace our ground truth observations with our predictions.There's an easy way to do this in python using the `np.roll(a, shift)` function.  That shifts everything in the array down by `shift`.  The function is **circular** so the value in element 0 is moved to be the final value in the array.<jupyter_code>y_pred = 999
current_X = np.array([1, 2, 3, 4])
current_X = np.roll(current_X, shift=-1)
current_X
current_X = np.array([1, 2, 3, 4])
current_X = np.roll(current_X, shift=-1)
current_X[-1] = y_pred
current_X

def autoregressive_iterative_forecast(model, exog, h):
    '''
    h-step forecast for an autoregressive 
    model using the iterative prediction method.
    
    Conduct h one-step forecasts gradually
    replacing ground truth autoregressive X 
    values with predictions.
    
    Parameters:
    ------
    model: forecast object
        model that has a .predict(h) interface
        
    exog: array-like
        initial vector of lagged values (X)
    
    h: int
        forecast horizon. assumed to be > 0
    
    Returns:
    ------
    numpy.ndarray
        y_predictions
    '''
    y_preds = []
    current_X = exog
    for i in range(h):
        y_pred = model.predict(current_X)[0]
        y_preds.append(y_pred)

        current_X = np.roll(current_X, shift=-1)
        current_X[-1] = y_pred

    return np.array(y_preds)
H = 5
y_preds = autoregressive_iterative_forecast(model, X_test[0], h=H)
print(f'Iterative forecast: {y_preds}')
print(f'Ground truth y: {y_test[:H].T[0]}')<jupyter_output>Iterative forecast: [-0.66561346 -0.50407924 -0.32244898 -0.12796369  0.07162311]
Ground truth y: [-0.66561346 -0.50407924 -0.32244898 -0.12796369  0.07162311]<jupyter_text>Adding some noise To make this a bit more interesting we will add some normally distributed noise to the synthetic time series.<jupyter_code>#set the random seed so that we all get the same results
np.random.seed(12)
t = np.arange(200)
ts_data = np.cos(0.2 * t)
noise = np.random.normal(loc=0.0, scale=0.3, size=200)
ts_data = ts_data + noise
plt.plot(ts_data);
def select_model_via_aic(train, test, window_sizes):  
    best = {'aic':np.Inf,
            'model': None,
            'ws': None}
    
    for ws in window_sizes:
    
        #preprocess time series training and test sets
        X_train, y_train = sliding_window(train, window_size=ws)
        X_test, y_test = sliding_window(test, window_size=ws)

        X_train = sm.add_constant(X_train)
        model = OLS(endog=y_train, exog=X_train).fit()
        
        if model.aic < best['aic']:
            best['aic'] = model.aic
            best['model'] = model
            best['ws'] = ws
            
    return best
ws = [i for i in range(2, 20, 2)]
print(f'Window sizes to test {ws}')

#train test split 150 and 50
train, test = ts_data[:150], ts_data[150:]

best = select_model_via_aic(train, test, ws)
print(f"lags included: {best['ws']}")
plt.plot(ts_data[best['ws']+1:], label='ground truth')
plt.plot(best['model'].fittedvalues, label='AR OLS fitted')
plt.legend();
#setup test data
ws = best['ws']

#preprocess time series training and test sets
X_train, y_train = sliding_window(train, window_size=ws)
X_test, y_test = sliding_window(test, window_size=ws)

X_train = sm.add_constant(X_train)
X_test = sm.add_constant(X_test)

model = OLS(endog=y_train, exog=X_train).fit()

#make iterative predictions
H = len(y_test)
y_preds_iter = autoregressive_iterative_forecast(model, X_test[0], h=H)

#plot
plt.plot(y_preds_iter, label='iterative forecast method')
plt.plot(y_test, label='ground truth')
plt.legend()<jupyter_output><empty_output><jupyter_text>The direct h-step forecasting method.In the direct method to forecast h-steps ahead we have **$h$ forecasting models**.  Each model provides a single point forecast from a step ahead.  In the example here, y_test is of length 31 periods.  The direct method requires 31 OLS models to make its prediction!Recall the `sliding_window` function.  We ignored an optional parameter `horizon` in the iterative example.  By default `horizon=1` i.e. the function returns target values that are only a single period ahead.  We can vary the step size by increasing the value of horizon.  Let's start by assigning `horizon` the value 2.<jupyter_code># creating the X and y training data
X_train, y_train = sliding_window(train, window_size=2, horizon=2)
#print the first 5 values in the original time series.  
#Note that y_train[0] is 2 steps ahead of X_train[0]
print(f'ts_data[:5]: {ts_data[:5]}')
print(f'X_train[0]: {X_train[0]}')
print(f'y_train[0]: {y_train[0]}')<jupyter_output>ts_data[:5]: [1.14189575 0.77563881 0.99379284 0.31511492 0.92264956]
X_train[0]: [1.14189575 0.77563881]
y_train[0]: [0.31511492]<jupyter_text>If we assign `horizon` the value 3 we get the same behaviour, but `y_train[0]` is now 3 steps ahead of `X_train[0]`<jupyter_code>X_train, y_train = sliding_window(train, window_size=2, horizon=3)
print(f'ts_data[:5]: {ts_data[:5]}')
print(f'X_train[0]: {X_train[0]}')
print(f'y_train[0]: {y_train[0]}')<jupyter_output>ts_data[:5]: [1.14189575 0.77563881 0.99379284 0.31511492 0.92264956]
X_train[0]: [1.14189575 0.77563881]
y_train[0]: [0.92264956]<jupyter_text>**Training multiple models**For OLS it is simple and quick to train the 31 models.  1. Create a for loop and set it to iterate 31 times. 2. In each loop call `sliding_window` setting `horizon` to the iteration number + 13. Train the model and save in a list.<jupyter_code>models = []
horizon = 31

for h in range(horizon):
    X_train, y_train = sliding_window(train, window_size=18, horizon=h+1)
    X_train = sm.add_constant(X_train)
    model_h = OLS(endog=y_train, exog=X_train).fit() 
    models.append(model_h)<jupyter_output><empty_output><jupyter_text>We now create the `direct_forecast` function.  This is just a for loop to call the `.predict(exog)` method of each model.  Remember that the input to each model is **same** i.e. exog which in our case will be `X_test[0]`<jupyter_code>def direct_forecast(models, exog):
    '''
    h-step forecast for an autoregressive 
    model using the direct prediction method.
    
    Each model contained in @models has been trained
    to predict a unique number of steps ahead. 
    Each model forecasts and the results are 
    combined in an ordered array and returned.
    
    Parameters:
    ------
    models: list
        direct models each has has a .predict(exog) 
        interface
        
    exog: array-like
        initial vector of lagged values (X)
        
    Returns:
    ------
    numpy.ndarray
        y_predictions
    '''
    preds = []
    for model_h in models:
        pred_h = model_h.predict(exog=exog)
        preds.append(pred_h)
    
    return np.array(preds)
#for clarity let's recreate the X and y test data
X_test, y_test = sliding_window(test, window_size=18)
X_test = sm.add_constant(X_test)

#make the direct forecast
y_preds_direct = direct_forecast(models, X_test[0])
#plot the direct forecast against the test data
plt.plot(y_preds_direct, label='direct forecast method')
plt.plot(y_test, label='ground truth')
plt.legend()<jupyter_output><empty_output><jupyter_text>Like the iterative method the direct method looks a close match to the ground truth test set!  Let's plot all three datasets on the same chart and then take a look at the **RMSE** of each method.<jupyter_code>#plot iterative and direct
plt.plot(y_preds_direct, label='direct forecast method')
plt.plot(y_preds_iter, label='iterative forecast method')
plt.plot(y_test, label='ground truth')
plt.legend()
from statsmodels.tools.eval_measures import rmse
rmse(y_test, y_preds_iter.reshape(-1, 1))[0]
rmse(y_test, y_preds_direct)[0]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1032258064516129, 'fraction_numerical': 0.038117182356813696, 'mean_word_length': 4.332046332046332, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '636607', 'license': 'MIT', 'n_tokens_mistral': 5382, 'n_tokens_neox': 4993, 'n_words': 1686, 'path': '05_autoregression/autoregression1.ipynb', 'repo_name': 'health-data-science-OR/forecasting', 'size': '15190'}