# Prediction of a time series AR(1) vs AR(1) with exogenous variables vs Random Forest, why is the performance so different?

Extension of the previous question that only compare AR(1) vs. AR(1) with exogenous variables:

I am currently working on forecasting a time series y using three models: an AR(1) model and an AR(1) model with exogenous variables, and a Random Forest. The time series y is provided in the code, while the exogenous variables can be found in this Excel sheet: https://github.com/YoussufAbdelatif/cpb_forecast.

In my implementation, using a 70/30 split for training and testing, the AR(1) model and Random Forest perform reasonably well. However, the AR(1) model with exogenous variables shows a much better performance, with a significantly lower RMSE. The difference in RMSE is so pronounced that I am concerned there might be a coding issue. Is it possible that the inclusion of the exogenous variables alone can account for such a large difference, or could there be a mistake in my implementation? If yes, why cannot the Random Forest exploit that, as it can access the same information that the AR(1) with exognous variables uses.

Differences in RMSE:

RMSE of AR(1): 0.040171
RMSE of AR(1) with exogenous variables: 0.008977
RMSE of Random Forest: 0.034226


Plot of AR1 predictions vs. AR1 with exogenous values vs. Random Forest predictions vs. actual values:

Comment to the explanatory variables: The explanatory variables are already lagged, so for y in t=0, X is in t-1 + The explanatory variables contain four lags of the targeted variable.

Code:

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy
import sklearn as sk
import plotly.graph_objects as go
import tensorflow as tf
import torch
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.neighbors import KNeighborsClassifier
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import re
from datetime import datetime, timedelta
from dateutil.relativedelta import relativedelta
from sklearn.linear_model import ElasticNet, ElasticNetCV, enet_path, lasso_path, LinearRegression
from sklearn.datasets import make_regression
from sklearn.model_selection import TimeSeriesSplit, RandomizedSearchCV, train_test_split
from itertools import product, cycle
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, LSTM
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import FormatStrFormatter
import xgboost as xgb
from xgboost import XGBClassifier, XGBRegressor
import statsmodels.api as sm
from statsmodels.formula.api import ols
from skglm import GeneralizedLinearEstimator
from tensorflow.keras.layers import LSTM, Dropout, Dense
from kerastuner.tuners import BayesianOptimization
import shutil
from statsmodels.tsa.api import VAR
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.arima.model import ARIMA

# Define the time series data
y = [
0.0250477, 0.0251129, 0.0356179, 0.0307251, 0.0308406, 0.0261657, 0.0202818,
-0.00931811, -0.0220432, -0.0189889, -0.00261438, 0.0192056, 0.0235917, 0.0205837,
0.00866395, 0.0147383, 0.00198899, 0.0171268, 0.0441156, 0.0226473, 0.0311145,
0.00671943, 0.0248557, 0.0060732, 0.030977, 0.0050077, 0.0341127, 0.0267606,
0.0180492, 0.00875824, 0.0247812, 0.0140976, 0.00855742, 0.0136495, 0.01244,
0.0153917, -0.00566426, -0.00313956, -0.0786389, -0.107904, 0.00624136, 0.0472657,
0.0483937, 0.0304634, 0.0378904, 0.0125033, 0.0229479, 0.0131371, -0.00232573,
0.0105381, -0.00420288, 0.00491876, 0.0110841, 0.00181148, -0.00317995, 0.0130418,
0.000957056, 0.00511998, 0.0081932, 0.00593517, 0.00472012, 0.0121062, 0.01092,
0.00214864, -0.00760691, 0.00852348, 0.00707222, -0.00477883, 0.00401124, 0.00676562,
0.0190308, 0.0138147, 0.00975719, 0.0101072, 0.0170982, 0.00335135, 0.00595302,
0.00953274, -0.00623326, -0.00880799, -0.000888626, 0.00423146, -0.00450886, -0.0395503,
-0.105516, 0.120913, 0.0439542, 0.0209713, 0.0111634, -0.00319514, 0.0330441,
-0.00239651, 0.00381834, 0.0112099, -0.0192566, -0.00970307, -0.00110293, -0.00331245,
0.00270514
]

# Load exogenous variables from an Excel file

# Convert y from list to numpy array with dtype=object
y = np.array(y, dtype=object)

# Convert y to numeric and ensure it is of type float64
y = pd.to_numeric(y, errors='coerce').astype(np.float64)

# Convert X to numeric types, handling errors and ensuring all data is float64
X = X.apply(pd.to_numeric, errors='coerce')
X = X.to_numpy()

# Split the data into training and testing sets (70/30 split)
split_index = int(len(X) * 0.7)
X_train = X[:split_index]    # Training features
X_test = pd.DataFrame(X[split_index:])  # Testing features (converted back to DataFrame for easier handling)
y_train = y[:split_index]    # Training target variable
y_test = y[split_index:]     # Testing target variable

# Initialize arrays for predictions
lr_preds = np.zeros(len(y_test))  # Predictions for Linear Regression with exogenous variables
ar1_preds = np.zeros(len(y_test))  # Predictions for AR(1) model
rf_preds = np.zeros(len(y_test))

rf_model = RandomForestRegressor(random_state=42,n_estimators = 1000,n_jobs = -1)

for i in range(len(y_test)):
print(i, (len(X_train) + i), len(y_test))

# Prepare training data for this iteration
x_train_data = X[:len(X_train) + i, :]  # Features up to the current point
y_train_data = y[:len(y_train) + i]     # Target values up to the current point

# Estimate AR(1) model (AutoReg model in statsmodels)
if len(y_train_data) > 1:
ar1_model = ARIMA(y_train_data, order=(1, 0, 0))  # AR(1) model (AutoReg model with ARIMA)
ar1_fit = ar1_model.fit()                       # Fit the model
ar1_preds[i] = ar1_fit.forecast()[0]           # Forecast the next value

# Estimate Linear Regression model with exogenous variables
if len(y_train_data) > 1:
lr_model = AutoReg(y_train_data, lags=1, exog=x_train_data, old_names=False).fit()  # Fit AutoReg model with exogenous variables
lr_preds[i] = lr_model.predict(start=len(y_train_data), end=len(y_train_data), exog_oos=np.array(X_test.iloc[i, :]).reshape(1, -1))[0]  # Predict the next value

# Random Forest prediction
if len(y_train_data) > 1:
rf_model.fit(x_train_data, y_train_data.ravel())
forecast_rf = rf_model.predict(np.array(X_test.iloc[i, :]).reshape(1, -1))
rf_preds[i] = forecast_rf

# Calculate RMSE for both models
rmse_ar1 = np.sqrt(mean_squared_error(y_test, ar1_preds))  # RMSE for AR(1) model
rmse_lr = np.sqrt(mean_squared_error(y_test, lr_preds))    # RMSE for AR(1) with exogenous variables
rmse_rf = np.sqrt(mean_squared_error(y_test, rf_preds))

# Print RMSE values
print(f"RMSE of AR(1): {rmse_ar1:.6f}")
print(f"RMSE of AR(1) with exogenous variables: {rmse_lr:.6f}")
print(f"RMSE of Random Forest: {rmse_rf:.6f}")

# Plot actual vs predicted values
plt.figure(figsize=(14, 7))
plt.plot(y_test, label='Actual')  # Plot actual values
plt.plot(lr_preds, label='AR(1) with exogenous variables')  # Plot predictions from AR(1) with exogenous variables
plt.plot(ar1_preds, label='AR(1) Predictions')  # Plot predictions from AR(1) model
plt.plot(rf_preds, label='Random Forest Predictions')  # Plot predictions from AR(1) model
plt.title('Actual vs Predicted Values')  # Plot title
plt.legend()  # Show legend
plt.show()  # Display the plot

• Whatever your data is, it looks like your exogeneous features are doing a great job picking up, explaining and forecasting the target. What are you concerned about? (Of course, in a "production" setting, you would need to forecast the features themselves. But this does not seem like what you are worrying about at this point.) Commented Jul 31 at 10:15
• @StephanKolassa The features I'm using are already lagged, so I'm not concerned about forecasting them. What I find puzzling is that, despite fitting a Random Forest model to the data, it does not capture the variability as well as the ARMAX model. This is surprising since the Random Forest also uses the same explanatory variables and lags of y. I think that I will add this to the question. Commented Jul 31 at 11:23
• Ah. I read your question as comparing an AR(1) to an AR(1) with explanatory variables model. It would be good if you could focus your question (and your code) to the essentials. Commented Jul 31 at 11:25
• @StephanKolassa yes, this was the main point first. I thought that there might be somehow an issue of data leakage with the AR(1) with explanatory variables. But if this is not the case, my questoin extends to the relatively poorer performance of a random forest in this case. Commented Jul 31 at 11:32
• @StephanKolassa I have just edited the question. Commented Jul 31 at 11:41