Time series no train/test split

docs/_quarto.yml

@@ -289,9 +289,9 @@ website:
                     href: notes/predictive-modeling/regression/seasonality.qmd
                     text: "Seasonality Analysis"
 
-                  - section:
-                    href: notes/predictive-modeling/regression/autoregressive-models.qmd
-                    text: "Autoregressive Models"
+              - section:
+                href: notes/predictive-modeling/regression/autoregressive-models.qmd
+                text: "Autoregressive Models"
 
               - section:
                 href: notes/predictive-modeling/classification/index.qmd

docs/notes/predictive-modeling/regression/seasonality.qmd

@@ -1 +1,336 @@
 # Regression for Seasonality Analysis
+
+Let's explore an example of how to use regression to identify cyclical patterns and perform seasonality analysis with time series data.
+
+
+```{python}
+#| echo: false
+
+#import warnings
+#warnings.simplefilter(action='ignore', category=FutureWarning)
+
+from pandas import set_option
+set_option('display.max_rows', 6)
+```
+
+
+## Data Loading
+
+
+As an example time series dataset, let's consider this dataset of U.S. unemployment rates over time, from the Federal Reserve Economic Data (FRED).
+
+Fetching the data, going back as far as possible:
+
+```{python}
+from pandas_datareader import get_data_fred
+from datetime import datetime
+
+DATASET_NAME = "PAYNSA"
+df = get_data_fred(DATASET_NAME, start=datetime(1900,1,1))
+print(len(df))
+df
+```
+
+:::{.callout-tip title="Data Source"}
+Here is some more information about the ["PAYNSA" dataset](https://fred.stlouisfed.org/series/PAYNSA):
+
+> All Employees: Total Nonfarm, commonly known as Total Nonfarm Payroll, is a measure of the number of U.S. workers in the economy that excludes proprietors, private household employees, unpaid volunteers, farm employees, and the unincorporated self-employed.
+> This measure accounts for approximately 80 percent of the workers who contribute to Gross Domestic Product (GDP).
+> This measure provides useful insights into the current economic situation because it can represent the number of jobs added or lost in an economy.
+> Increases in employment might indicate that businesses are hiring which might also suggest that businesses are growing. Additionally, those who are newly employed have increased their personal incomes, which means (all else constant) their disposable incomes have also increased, thus fostering further economic expansion.
+> Generally, the U.S. labor force and levels of employment and unemployment are subject to fluctuations due to seasonal changes in weather, major holidays, and the opening and closing of schools.
+> The Bureau of Labor Statistics (BLS) adjusts the data to offset the seasonal effects to show non-seasonal changes: for example, women's participation in the labor force; or a general decline in the number of employees, a possible indication of a downturn in the economy. To closely examine seasonal and non-seasonal changes, the BLS releases two monthly statistical measures: the seasonally adjusted All Employees: Total Nonfarm (PAYEMS) and All Employees: Total Nonfarm (PAYNSA), which is not seasonally adjusted.
+:::
+
+
+
+
+Wrangling the data, including renaming columns and converting the date index to be datetime-aware, may make it easier for us to work with this data:
+
+```{python}
+from pandas import to_datetime
+
+df.rename(columns={DATASET_NAME: "employment"}, inplace=True)
+df.index.name = "date"
+df.index = to_datetime(df.index)
+df
+```
+
+
+
+## Data Exploration
+
+
+Visualizing the data:
+
+```{python}
+import plotly.express as px
+
+px.line(df, y="employment", height=450,
+        title="US Employment by month (non-seasonally adjusted)",
+        labels={"employment": "Employment (in thousands of persons)"},
+)
+```
+
+**Cyclical Patterns**
+
+Exploring cyclical patterns in the data:
+
+```{python}
+px.line(df[(df.index.year >= 1970) & (df.index.year <= 1980)], y="employment",
+        title="US Employment by month (selected years)", height=450,
+        labels={"Employment": "Employment (in thousands)"},
+)
+```
+
+:::{.callout-tip title="Interactive dataviz"}
+Hover over the dataviz to see which month(s) typically have higher employment, and which month(s) typically have lower employment.
+:::
+
+**Trend Analysis**
+
+Exploring trends:
+
+```{python}
+import plotly.express as px
+
+px.scatter(df, y="employment",  height=450,
+        title="US Employment by month (vs Trend)",
+        labels={"employment": "Employment (in thousands)"},
+        trendline="ols", trendline_color_override="red"
+)
+```
+
+Looks like evidence of a possible linear relationship. Let's perform a more formal regression analysis.
+
+
+## Data Encoding
+
+Because we need numeric features to perform a regression, we convert the dates to a linear time step of integers (after sorting the data first for good measure):
+
+```{python}
+df.sort_values(by="date", ascending=True, inplace=True)
+
+df["time_step"] = range(1, len(df) + 1)
+df
+```
+
+We will use the numeric time step as our input variable (`x`), to predict the employment (`y`).
+
+
+
+
+
+
+
+
+
+
+## Data Splitting
+
+**X/Y Split**
+
+Identifying dependent and independent variables:
+
+```{python}
+x = df[["time_step"]]
+
+y = df["employment"]
+
+print("X:", x.shape)
+print("Y:", y.shape)
+```
+
+**Adding Constants**
+
+We are going to use `statsmodels`, so we add a column of constant ones representing the intercept:
+
+```{python}
+import statsmodels.api as sm
+
+# adding in a column of constants, as per the OLS docs
+x = sm.add_constant(x)
+x.head()
+```
+
+**Train/Test Split**
+
+Splitting into training vs testing datasets:
+
+```{python}
+#from sklearn.model_selection import train_test_split
+#
+#x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=99)
+#print("TRAIN:", x_train.shape, y_train.shape)
+#print("TEST:", x_test.shape, y_test.shape)
+```
+
+Splitting data sequentially where earlier data is used in training and recent data is use for testing:
+
+```{python}
+#print(len(df))
+#
+#training_size = round(len(df) * .8)
+#print(training_size)
+#
+#x_train = x.iloc[:training_size] # slice all before
+#y_train = y.iloc[:training_size] # slice all before
+#
+#x_test = x.iloc[training_size:] # slice all after
+#y_test = y.iloc[training_size:] # slice all after
+#print("TRAIN:", x_train.shape)
+#print("TEST:", x_test.shape)
+```
+
+
+For this example, we will not split the data. To help illustrate a story about predictions over the entire time period.
+
+
+
+
+## Model Selection and Training
+
+Training a linear regression model on the training data:
+
+```{python}
+import statsmodels.api as sm
+
+#model = sm.OLS(y_train, x_train, missing="drop")
+model = sm.OLS(y, x, missing="drop")
+print(type(model))
+
+results = model.fit()
+print(type(results))
+```
+
+Examining training results:
+
+```{python}
+print(results.summary())
+```
+
+
+```{python}
+print(results.params)
+print("------------")
+print(f"y = {results.params['time_step'].round(3)}x + {results.params['const'].round(3)}")
+```
+
+```{python}
+df["prediction"] = results.fittedvalues
+df["residual"] = results.resid
+```
+
+```{python}
+#from pandas import DataFrame
+#
+## get all rows from the original dataset that wound up in the training set:
+#training_set = df.loc[x_train.index].copy()
+#print(len(training_set))
+#
+## create a dataset for the predictions and the residuals:
+#training_preds = DataFrame({
+#    "prediction": results.fittedvalues,
+#    "residual": results.resid
+#})
+## merge the training set with the results:
+#training_set = training_set.merge(training_preds,
+#    how="inner", left_index=True, right_index=True
+#)
+#
+## calculate error for each datapoint:
+#training_set
+```
+
+**Regression Trends**
+
+Plotting trend line:
+
+```{python}
+px.line(df, y=["employment", "prediction"], height=350,
+    title="US Employment (monthly) vs linear trend",
+    labels={"value":""}
+)
+```
+
+**Regression Residuals**
+
+Plotting residuals:
+
+
+```{python}
+px.line(df, y="residual",
+    title="US Employment (monthly) vs linear trend residuals", height=350
+)
+```
+
+#### Seasonality via Means of Periodic Residuals
+
+Observe there may be some cyclical patterns in the residuals, by calculating periodic means:
+
+```{python}
+#| echo: false
+
+from pandas import set_option
+set_option('display.max_rows', 15)
+```
+
+```{python}
+df = df.copy()
+df["year"] = df.index.year
+df["quarter"] = df.index.quarter
+df["month"] = df.index.month
+```
+
+```{python}
+df.groupby("quarter")["residual"].mean()
+```
+
+
+```{python}
+df.groupby("month")["residual"].mean()
+```
+
+```{python}
+#| echo: false
+
+from pandas import set_option
+set_option('display.max_rows', 6)
+```
+
+
+
+#### Seasonality via Regression on Periodic Residuals
+
+
+```{python}
+# https://pandas.pydata.org/docs/reference/api/pandas.get_dummies.html
+# "one hot encode" the monthly values:
+from pandas import get_dummies as one_hot_encode
+
+x_monthly = one_hot_encode(df["month"])
+x_monthly.columns=["Jan", "Feb", "Mar", "Apr",
+                "May", "Jun", "Jul", "Aug",
+                "Sep", "Oct", "Nov", "Dec"]
+x_monthly = x_monthly.astype(int)
+x_monthly
+```
+
+
+Predict the residual (i.e. degree to which we will be over or under trend), based on which month it is?
+
+```{python}
+
+
+y_monthly = df["residual"]
+
+ols_monthly = sm.OLS(y_monthly, x_monthly)
+print(type(ols_monthly))
+
+results_monthly = ols_monthly.fit()
+print(type(results_monthly))
+
+print(results_monthly.summary())
+```
+
+Observe the coefficients tell us how each month contributes towards the regression residuals, in other words, for each month, to what degree does the model predict we will be above or below trend?