Time Series Forecasting with Prophet, NeuralProphet, and LSTM

Table of Contents

Introduction

Time series forecasting is one of the most impactful applications of data science in industry. From predicting stock prices and energy consumption to forecasting product demand and server load, the ability to anticipate future values from historical patterns is invaluable.

This tutorial covers three powerful approaches to time series forecasting: Facebook Prophet for robust additive decomposition, NeuralProphet for neural network-enhanced forecasting with autoregression, and LSTM (Long Short-Term Memory) networks built in PyTorch for capturing complex sequential dependencies.

By the end of this tutorial, you will be able to build, evaluate, and compare production-grade forecasting models using all three frameworks.

Prerequisites

• Python 3.9+
• Basic understanding of statistics (mean, variance, correlation)
• Familiarity with pandas and matplotlib

Install the required packages:

pip install prophet neuralprophet torch pandas numpy matplotlib scikit-learn

Time Series Fundamentals

A time series is a sequence of data points indexed in time order. Before building any model, you must understand the key components that make up a time series signal.

Components of a Time Series

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

Generate a synthetic time series with known components
np.random.seed(42)
dates = pd.daterange(start="2020-01-01", periods=730, freq="D")


trend = np.linspace(50, 150, 730)
seasonality = 20  np.sin(2  np.pi  np.arange(730) / 365.25)

noise = np.random.normal(0, 5, 730)

y = trend + seasonality + noise

df = pd.DataFrame({"ds": dates, "y": y})

fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(dates, y, color="steelblue")
axes[0].settitle("Observed Time Series")
axes[1].plot(dates, trend, color="darkorange")
axes[1].settitle("Trend Component")

axes[2].plot(dates, seasonality, color="green")
axes[2].settitle("Seasonal Component")
axes[3].plot(dates, noise, color="red", alpha=0.6)
axes[3].settitle("Noise / Residual")

plt.tightlayout()
plt.savefig("timeseriesdecomposition.png", dpi=150)
plt.show()

from statsmodels.tsa.stattools import adfuller

result = adfuller(df["y"].values)
print(f"ADF Statistic: {result[0]:.4f}")
print(f"p-value: {result[1]:.4f}")

if result[1] < 0.05:
    print("The series is stationary (reject null hypothesis).")
else:
    print("The series is non-stationary. Consider differencing.")

splitdate = "2021-07-01"
train = df[df["ds"] < splitdate].copy()
test = df[df["ds"] >= splitdate].copy()


print(f"Training samples: {len(train)}")
print(f"Testing samples:  {len(test)}")

from prophet import Prophet

model = Prophet(
    growth="linear",
    changepointpriorscale=0.05,   # flexibility of trend changes

    seasonalitypriorscale=10.0,   # strength of seasonality

    yearlyseasonality=True,
    weeklyseasonality=False,

    dailyseasonality=False,
)
model.fit(train)

Create future dataframe for predictions
future = model.makefuturedataframe(periods=len(test), freq="D")
forecast = model.predict(future)

Visualize forecast
fig = model.plot(forecast)
plt.title("Prophet Forecast")
plt.savefig("prophetforecast.png", dpi=150)

plt.show()

holidays = pd.DataFrame({
    "holiday": "companyevent",
    "ds": pd.todatetime(["2020-12-25", "2021-12-25", "2020-01-01", "2021-01-01"]),

    "lowerwindow": -1,
    "upperwindow": 1,

})

modelwithholidays = Prophet(holidays=holidays)

modelwithholidays.fit(train)

forecastholidays = modelwithholidays.predict(future)

fig = modelwithholidays.plotcomponents(forecastholidays)
plt.savefig("prophetcomponents.png", dpi=150)

plt.show()

# Add monthly seasonality
modelcustom = Prophet()
modelcustom.addseasonality(name="monthly", period=30.5, fourierorder=5)


Add an external regressor (e.g., temperature)
train["temperature"] = np.random.normal(25, 5, len(train))
test["temperature"] = np.random.normal(25, 5, len(test))

modelcustom.addregressor("temperature")

modelcustom.fit(train)

futurecustom = pd.concat([train, test])[["ds", "temperature"]]

forecastcustom = modelcustom.predict(futurecustom)

from prophet.diagnostics import crossvalidation, performancemetrics

dfcv = crossvalidation(
    model,
    initial="365 days",
    period="90 days",
    horizon="90 days",
)

dfmetrics = performancemetrics(dfcv)

print(dfmetrics[["horizon", "mae", "rmse", "mape"]].tail(10))

from neuralprophet import NeuralProphet

npmodel = NeuralProphet(

    growth="linear",
    nchangepoints=10,
    yearlyseasonality=True,

    weeklyseasonality=False,
    dailyseasonality=False,

    nforecasts=1,
    nlags=0,              # no autoregression yet

    learningrate=0.1,
    epochs=100,
    batchsize=64,

)

metrics = npmodel.fit(train, freq="D")
forecastnp = npmodel.predict(test)

npmodel.plot(forecastnp)
plt.title("NeuralProphet Forecast")
plt.savefig("neuralprophetforecast.png", dpi=150)

plt.show()

nparmodel = NeuralProphet(
    growth="linear",
    nchangepoints=10,
    yearlyseasonality=True,

    nforecasts=30,        # predict 30 steps ahead
    nlags=60,             # use 60 past values as input

    learningrate=0.01,
    epochs=200,
    batchsize=64,

    arlayers=[32, 16],    # hidden layers for AR-Net
)

metricsar = nparmodel.fit(train, freq="D", validationdf=test)

forecastar = nparmodel.predict(test)

nparmodel.plot(forecastar)
plt.title("NeuralProphet AR Forecast (30-step)")
plt.savefig("neuralprophetarforecast.png", dpi=150)
plt.show()

train["temp"] = np.random.normal(25, 5, len(train))
test["temp"] = np.random.normal(25, 5, len(test))

npregmodel = NeuralProphet(
    nforecasts=7,

    nlags=30,
    yearlyseasonality=True,

    epochs=150,
)
npregmodel.addlaggedregressor("temp", nlags=14)

metricsreg = npregmodel.fit(train, freq="D", validationdf=test)
forecastreg = npregmodel.predict(test)

import torch
from torch import nn
from torch.utils.data import Dataset, DataLoader
from sklearn.preprocessing import MinMaxScaler

scaler = MinMaxScaler(featurerange=(0, 1))
scaledvalues = scaler.fittransform(df["y"].values.reshape(-1, 1))

def createsequences(data, seqlength):
    xs, ys = [], []
    for i in range(len(data) - seqlength):

        xs.append(data[i : i + seqlength])
        ys.append(data[i + seqlength])

    return np.array(xs), np.array(ys)

SEQLENGTH = 60
X, y = createsequences(scaledvalues, SEQLENGTH)


splitidx = int(len(X)  0.8)
Xtrain, Xtest = X[:splitidx], X[splitidx:]
ytrain, ytest = y[:splitidx], y[splitidx:]

Xtraint = torch.FloatTensor(Xtrain)

ytraint = torch.FloatTensor(ytrain)
Xtestt = torch.FloatTensor(Xtest)

ytestt = torch.FloatTensor(ytest)

class TimeSeriesDataset(Dataset):
    def init(self, X, y):
        self.X = X
        self.y = y

    def len(self):
        return len(self.X)

    def getitem(self, idx):
        return self.X[idx], self.y[idx]

traindataset = TimeSeriesDataset(Xtraint, ytraint)

trainloader = DataLoader(traindataset, batchsize=32, shuffle=True)

LSTM Model Architecture

class LSTMForecaster(nn.Module):
    def init(self, inputsize=1, hiddensize=64, numlayers=2, dropout=0.2):

        super(LSTMForecaster, self).init()
        self.hiddensize = hiddensize

        self.numlayers = numlayers


        self.lstm = nn.LSTM(
            inputsize=inputsize,

            hiddensize=hiddensize,

            numlayers=numlayers,

            batchfirst=True,
            dropout=dropout,
        )
        self.fc = nn.Sequential(
            nn.Linear(hiddensize, 32),

            nn.ReLU(),
            nn.Linear(32, 1),
        )

    def forward(self, x):
        # x shape: (batch, seqlength, inputsize)

        lstmout, (hn, cn) = self.lstm(x)
        # Use the last hidden state
        lasthidden = lstmout[:, -1, :]
        output = self.fc(lasthidden)

        return output

device = torch.device("cuda" if torch.cuda.isavailable() else "cpu")
model = LSTMForecaster(inputsize=1, hiddensize=64, numlayers=2, dropout=0.2).to(device)

print(model)

criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)
scheduler = torch.optim.lrscheduler.ReduceLROnPlateau(optimizer, patience=5, factor=0.5)

numepochs = 100

bestvalloss = float("inf")

trainlosses, vallosses = [], []


for epoch in range(numepochs):
    model.train()
    epochloss = 0.0

    for batchX, batchy in trainloader:
        batchX, batchy = batchX.to(device), batchy.to(device)
        optimizer.zerograd()

        predictions = model(batchX)
        loss = criterion(predictions, batchy)

        loss.backward()
        torch.nn.utils.clipgradnorm(model.parameters(), maxnorm=1.0)

        optimizer.step()
        epochloss += loss.item()

    avgtrainloss = epochloss / len(trainloader)
    trainlosses.append(avgtrainloss)


    # Validation
    model.eval()
    with torch.nograd():
        valpred = model(Xtestt.to(device))

        valloss = criterion(valpred, ytestt.to(device)).item()

        vallosses.append(valloss)


    scheduler.step(valloss)

    if valloss < bestvalloss:

        bestvalloss = valloss
        torch.save(model.statedict(), "bestlstmmodel.pt")


    if (epoch + 1) % 10 == 0:
        print(f"Epoch [{epoch+1}/{numepochs}] Train Loss: {avgtrainloss:.6f} | Val Loss: {valloss:.6f}")


Plot learning curves
plt.figure(figsize=(10, 5))
plt.plot(trainlosses, label="Train Loss")
plt.plot(vallosses, label="Validation Loss")

plt.xlabel("Epoch")
plt.ylabel("MSE Loss")
plt.legend()
plt.title("LSTM Training Curves")
plt.savefig("lstmtrainingcurves.png", dpi=150)

plt.show()

model.loadstatedict(torch.load("bestlstmmodel.pt"))
model.eval()

with torch.nograd():
    testpredictions = model(Xtestt.to(device)).cpu().numpy()


Inverse transform to original scale
testpredictionsinv = scaler.inversetransform(testpredictions)

ytestinv = scaler.inversetransform(ytest)


plt.figure(figsize=(14, 5))
plt.plot(ytestinv, label="Actual", color="steelblue")

plt.plot(testpredictionsinv, label="LSTM Prediction", color="darkorange")

plt.title("LSTM Forecast vs Actual")
plt.legend()
plt.savefig("lstmforecast.png", dpi=150)
plt.show()

Anomaly Detection in Time Series

Forecasting models can double as anomaly detectors. Points that deviate significantly from the forecast confidence interval are flagged as anomalies.

Prophet-Based Anomaly Detection

forecastfull = model.predict(model.makefuturedataframe(periods=0))
dfmerged = df.merge(forecastfull[["ds", "yhat", "yhatlower", "yhatupper"]], on="ds")


dfmerged["anomaly"] = (dfmerged["y"] < dfmerged["yhatlower"]) | (dfmerged["y"] > dfmerged["yhatupper"])

anomalies = dfmerged[dfmerged["anomaly"]]
print(f"Detected {len(anomalies)} anomalies out of {len(dfmerged)} data points.")


plt.figure(figsize=(14, 5))
plt.plot(dfmerged["ds"], dfmerged["y"], label="Actual", alpha=0.7)

plt.fillbetween(dfmerged["ds"], dfmerged["yhatlower"], dfmerged["yhatupper"], alpha=0.2, color="blue")

plt.scatter(anomalies["ds"], anomalies["y"], color="red", s=50, label="Anomaly", zorder=5)
plt.legend()
plt.title("Anomaly Detection with Prophet")
plt.savefig("anomalydetection.png", dpi=150)
plt.show()

Z-Score Residual Method

residuals = dfmerged["y"] - dfmerged["yhat"]
zscores = (residuals - residuals.mean()) / residuals.std()


threshold = 2.5
dfmerged["zanomaly"] = np.abs(zscores) > threshold
print(f"Z-score anomalies: {dfmerged['zanomaly'].sum()}")

Evaluation Metrics

Always evaluate your forecasting models with multiple metrics to get a complete picture.

from sklearn.metrics import meanabsoluteerror, meansquarederror

def evaluateforecast(actual, predicted, modelname="Model"):
    mae = meanabsoluteerror(actual, predicted)
    rmse = np.sqrt(meansquarederror(actual, predicted))
    # Avoid division by zero in MAPE
    mask = actual != 0
    mape = np.mean(np.abs((actual[mask] - predicted[mask]) / actual[mask]))  100


    print(f"--- {modelname} ---")

    print(f"  MAE:  {mae:.4f}")
    print(f"  RMSE: {rmse:.4f}")
    print(f"  MAPE: {mape:.2f}%")
    return {"mae": mae, "rmse": rmse, "mape": mape}

Example usage
actualvalues = test["y"].values
prophetpreds = forecast[forecast["ds"].isin(test["ds"])]["yhat"].values


metricsprophet = evaluateforecast(actualvalues, prophetpreds, "Prophet")

class MultiStepLSTM(nn.Module):
    def init(self, inputsize=1, hiddensize=64, numlayers=2, forecasthorizon=7):

        super(MultiStepLSTM, self).init()
        self.lstm = nn.LSTM(inputsize, hiddensize, numlayers, batchfirst=True, dropout=0.2)

        self.fc = nn.Linear(hiddensize, forecasthorizon)


    def forward(self, x):
        lstmout,  = self.lstm(x)

        lasthidden = lstmout[:, -1, :]

        return self.fc(lasthidden)

This model outputs 7 values at once
multimodel = MultiStepLSTM(forecasthorizon=7).to(device)

Recursive Multi-Step

def recursiveforecast(model, initialsequence, steps, scaler, device):
    model.eval()
    predictions = []
    currentseq = initialsequence.clone()

    for  in range(steps):

        with torch.nograd():
            pred = model(currentseq.unsqueeze(0).to(device))

        predictions.append(pred.cpu().item())
        # Shift the sequence and append the new prediction
        newval = torch.FloatTensor([[pred.cpu().item()]])
        currentseq = torch.cat([currentseq[1:], newval], dim=0)


    predictions = np.array(predictions).reshape(-1, 1)
    return scaler.inversetransform(predictions)

Forecast 30 steps ahead
initialseq = Xtestt[0]

longforecast = recursiveforecast(model, initialseq, steps=30, scaler=scaler, device=device)

Model Comparison

results = {}

Prophet
prophetpred = forecast[forecast["ds"].isin(test["ds"])]["yhat"].values

results["Prophet"] = evaluateforecast(test["y"].values, prophetpred, "Prophet")


NeuralProphet
nppred = forecastnp["yhat1"].values[:len(test)]

results["NeuralProphet"] = evaluateforecast(test["y"].values, nppred, "NeuralProphet")


LSTM
results["LSTM"] = evaluateforecast(ytestinv.flatten(), testpredictionsinv.flatten(), "LSTM")

Summary table
comparisondf = pd.DataFrame(results).T

comparisondf.columns = ["MAE", "RMSE", "MAPE"]
print("\n=== Model Comparison ===")
print(comparisondf.tostring())

Visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, metric in enumerate(["MAE", "RMSE", "MAPE"]):
    axes[i].bar(comparisondf.index, comparisondf[metric], color=["#4C72B0", "#55A868", "#C44E52"])
    axes[i].settitle(metric)

    axes[i].setylabel(metric)
plt.tightlayout()

plt.savefig("modelcomparison.png", dpi=150)
plt.show()

Best Practices

Conclusion

In this tutorial, we explored three complementary approaches to time series forecasting:

• Prophet excels at interpretable, decomposition-based forecasting with minimal tuning. It is ideal for business time series with clear trend and seasonal patterns.
• NeuralProphet extends Prophet with autoregressive neural components, making it suitable for series where recent past values are strong predictors.
• LSTM provides maximum flexibility for complex, nonlinear patterns but requires more data preparation and hyperparameter tuning.

The best approach depends on your data characteristics, accuracy requirements, and operational constraints. For production systems, consider starting with Prophet for a strong baseline, then incrementally adding complexity with NeuralProphet or LSTM if needed. Always validate with proper time-based cross-validation and multiple evaluation metrics.