HFW 15: Trend Interlude

We switch to trend forecasting models which perform almost as well as the benchmarks and easily explainable too

Today is as good a day for a transition as any other. Over the last 15 posts (we used zero-based indexing since we're soo Pythonic!), we've examined various univariate-based models—from naive to autoregressive to ARIMA to GARCH—to forecast the revenues for 50 companies representing ten S&P 500 sectors. Performance varied meaningfully, but the naïve and autoregressive growth rate (AR%(1)) models topped the list. As noted in our second full post, we're using the values of the variable itself— lagged, log-shifted, or residualed (our own neologism)—to forecast the future. While this has worked to with varying degrees of success, the point was we would be using endogenous variables, rather than exogenous ones to test the different algorithms. That means, we're not supplying variables like interest rates, GDP, or Google searches to the model.

What about time? Is it an endogenous or exogenous variable? Pedantically, it is exogenous because it's fixed and (theoretically) outside of the system being modeled. Commonsensically, it's endogenous. Time is not fixed. Time-series data wouldn't be such were it not for time. (As a side comment, we would argue all data are time series data, but that is a philosophical argument we'll save for another time.) If your model to forecast time-series data uses only the data to be forecast, then adding time as a variable is indeed exogenous. But time is implicit in all time-series forecasts. So why not forecast using time as an explicit variable? The result: trend forecasting!

For fundamental equity analysts this may seem like the most convoluted path to something that is pretty intuitive—trendline your forecast. But we had to crawl before we could walk. For the next few posts, we'll be using the sktime package. Anyone who has used scikit-learn (hard to imagine a data scientist that hasn't!) will be familiar with the .fit() and .predict() methods used in sktime. This makes it easy to get up and running with a ton of algorithms at one's fingertips along with many of pipelines that makes scikit-learn user-friendly and production-ready.

Ok, this isn't an infomercial. Let's get on with the data. We set up a basic trend forecast for all the revenue data we have. That is, we build the model using 16 quarters of data, with the independent variable being an integer time step. We then forecast the next four quarters as the next four sequential time steps. This is essentially a linear regression with integer time steps as the x variables and the quarterly revenues as the y variables. Let's show the scaled root mean-squared error chart that is now so common.

Those who have been reading each post with undivided attention will note immediately how well the trend model performs on companies in XLK, the technology ETF. This, of course, includes the darling of AI, Nvidia, whose meteoric revenue growth is nicely captured by a trend model.

Let's see how it performs relative to benchmarks.

Definitely, an outperformer on XLK. But besides that, XLP and XLC, the trend model is often worse than the naïve benchmark. Nonetheless, differences are modest, in aggregate, as we show below.

Indeed, the differences in scaled RMSE from best (naive) to worst (AR(1)) are less than 2% points.

What's the takeaway? Trend models yield pretty good forecasts that are eminently explainable. If you need to explain your forecasts to non-technical stakeholders, which would you prefer using? A trend model or a five-minute detour into the foundations of autoregressive integrated moving averages? Let us know in the comments!

Code below.

# Load packages
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
import statsmodels.api as sm
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import matplotlib as mpl
import pickle
from sktime.forecasting.trend import TrendForecaster
from sktime.utils.plotting import plot_series

# Assign chart style
plt.style.use('seaborn-v0_8')
plt.rcParams["figure.figsize"] = (12,6)

# Handy functions
# Functions to load and save files
def save_dict_to_file(data, filename):
    with open(filename, 'wb') as f:
        pickle.dump(data, f)

def load_dict_from_file(filename):
    with open(filename, 'rb') as f:
        return pickle.load(f)

# Create functions for indexing
def create_index(series):
    if series.iloc[0] > 0:
        return series/series.iloc[0] * 100
    else:
        return (series - series.iloc[0])/-series.iloc[0] * 100

# Function for scaled RMSE
def get_rmse_scaled(series):
    return np.sqrt(np.mean((series['actual']-series['predicted'])**2))/np.mean(series['actual'])

# Function to flatten
def flatten_df(dataf: pd.DataFrame, group_name:str, cols:list) -> pd.DataFrame:
    df_grouped = dataf.groupby(group_name)[cols].agg(list)

    for col in cols:
        df_grouped[col] = df_grouped[col].apply(lambda x: np.concatenate(x))
    
    df_long = df_grouped.apply(pd.Series.explode).reset_index()

    return df_long

# Symbols used
etf_symbols = ['XLF', 'XLI', 'XLE', 'XLK', 'XLV', 'XLY', 'XLP', 'XLB', 'XLU', 'XLC']

ticker_list = ["SHW", "LIN", "ECL", "FCX", "VMC",
               "XOM", "CVX", "COP", "WMB", "SLB",
               "JPM", "V", "MA", "BAC", "GS",
               "CAT", "RTX", "DE", "UNP", "BA",
               "AAPL", "MSFT", "NVDA", "ORCL", "CRM",
               "COST", "WMT", "PG", "KO", "PEP",
               "NEE", "D", "DUK", "VST", "SRE",
               "LLY", "UNH", "JNJ", "PFE", "MRK",
               "AMZN", "SBUX", "HD", "BKNG", "MCD",
               "META", "GOOG", "NFLX", "T", "DIS"
]

xlb = ["SHW", "LIN", "ECL", "FCX", "VMC"]
xle = ["XOM", "CVX", "COP", "WMB", "SLB"]
xlf = ["JPM", "V", "MA", "BAC", "GS"]
xli = ["CAT", "RTX", "DE", "UNP", "BA"]
xlk = ["AAPL", "MSFT", "NVDA", "ORCL", "CRM"]
xlp = ["COST", "WMT", "PG", "KO", "PEP"]
xlu = ["NEE", "D", "DUK", "VST", "SRE"]
xlv = ["LLY", "UNH", "JNJ", "PFE", "MRK"]
xly = ["AMZN", "SBUX", "HD", "BKNG", "MCD"]
xlc = ["META", "GOOG", "NFLX", "T", "DIS"]

sectors = [xlf,  xli,  xle,  xlk, xlv,  xly, xlp, xlb, xlu, xlc]

# Sector dictionary
sector_dict = {symbol.lower(): tickers for symbol, tickers in zip(etf_symbols, sectors)}
sector_map = {ticker.lower(): symbol for symbol, tickers in zip(etf_symbols, sectors) for ticker in tickers}

# Load data from disk
# See Code Walk-Throughs for how we built the data set
df_sector_dict = load_dict_from_file("path/to/data/simfin_df_rev_dict.pkl")

# Clean dataframes
df_rev_index_dict = {}
for key in sector_dict:
    temp_df = df_sector_dict[key].copy()
    col_1 = temp_df.columns[0]
    temp_df = temp_df[[col_1] + [x for x in temp_df.columns if 'revenue' in x]]
    temp_df.columns = ['date'] + [x.replace('revenue_', '').lower() for x in temp_df.columns[1:]]
    temp_idx = temp_df.copy()
    temp_idx[[x for x in temp_idx if 'date' not in x]] = temp_idx[[x for x in temp_idx if 'date' not in x]].apply(create_index)
    df_rev_index_dict[key] = temp_idx

# Create train/test dataframes
df_rev_train_dict = {}
df_rev_test_dict = {}
for key in df_rev_index_dict:
    # Create ticker list
    tickers = [x.lower() for x in sector_dict[key]]
    
    # Create dataframe fo all
    df_out = df_rev_index_dict[key]

    # Base train df
    df_train = df_out.loc[df_out['date'] < "2023-01-01"].copy()
    df_rev_train_dict[key] = df_train
    
    # Test df
    df_rev_test_dict[key] = df_out.loc[df_out['date'] >= "2023-01-01"]

# Create trend models
trend_models = {}
for key in df_rev_chg_train_dict:
    tickers = [x.lower() for x in sector_dict[key]]
    df_train = df_rev_train_dict[key]
    trend_models[key] = {}
    for ticker in tickers:
        try: 
            model = TrendForecaster().fit(df_train[ticker])
        except ValueError as e:
            print(f"{ticker} fails due to {e}")

        trend_models[key][ticker] = model

trend_err_df = pd.DataFrame(columns = ['sector', 'ticker', 'actual', 'predicted'])
count = 0
forecast_horizon = np.arange(1,5)
for key in trend_models:
    for ticker in trend_models[key]:
        # Forecast the next 4 quarters' volatility
        y_pred = trend_models[key][ticker].predict(forecast_horizon)
        y_act = df_rev_test_dict[key][ticker].values

        trend_err_df.loc[count] = [key.upper(), ticker, y_act, y_pred]
        count += 1

for ticker in trend_err_df['ticker'].to_list():
    actual = trend_err_df.loc[trend_err_df['ticker']==ticker, 'actual'].values[0]
    predicted = trend_err_df.loc[trend_err_df['ticker']==ticker, 'predicted'].values[0]
    rmse = np.sqrt(np.mean((actual - predicted)**2))
    rmse_scaled = rmse/actual.mean()
    trend_err_df.loc[trend_err_df['ticker']==ticker, 'rmse'] = rmse
    trend_err_df.loc[trend_err_df['ticker']==ticker, 'rmse_scaled'] = rmse_scaled
  
# Clean and flatten ar_err_df
trend_df_long = flatten_df(trend_err_df, 'sector', ['actual', 'predicted'])

# Graph scaled RMSE by sector
(trend_df_long.groupby('sector')[['actual', 'predicted']].apply(get_rmse_scaled).sort_values(ascending=True)*100).plot(kind='bar', rot=0)
plt.xlabel('')
plt.ylabel('RMSE (%)')
plt.title('Scaled RMSE by Sector for Trend Models')
plt.show()

comp_df = pd.read_pickle("hello_world/data/comp_new_df.pkl")
trend = trend_df_long.groupby('sector')[['actual', 'predicted']].apply(get_rmse_scaled).reset_index().rename(columns={0:'trend'})
comp_df = comp_df.merge(trend, how="left", on='sector')
comp_df.iloc[:, 1:].mean()*100

# Plot RMSE by sector and model
(comp_df.set_index('sector')[['naive', 'lr', 'chg','trend']].sort_values('naive')*100).plot(kind='bar', rot=0)
plt.xlabel('')
plt.ylabel('RMSE (%)')
plt.ylim(0, 110)
plt.legend(['Naive', 'AR(1)', 'AR%(1)', 'Trend'], loc='upper center', ncol=4)
plt.title('Scaled RMSE Comparisons by Sector and Model')
plt.show()

# Plot Average performance
x_lab_dict = {'naive': 'Naive', 'lr': 'AR(1)', 'chg': 'AR%(1)', 'trend':'Trend'}
x_labs = comp_df.iloc[:, 1:].mean().sort_values().index.to_list()

(comp_df.iloc[:, 1:].mean().sort_values()*100).plot(kind='bar', rot=0)
plt.xticks(np.arange(4), [x_lab_dict[x] for x in x_labs])
plt.ylabel('RMSE (%)')
plt.ylim(20,26)
plt.title('Average scaled RMSE by model')
plt.show()