# LSTM vf updated lstm only clean v5 populated-1

course: Module 3 — Deep Learning & NLP
module: Module-3-Deep-Learning-NLP
type: notebook
source_url: https://personal-learn.armco.dev/files/Module-3-Deep-Learning-NLP/General/Lab_Materials-28-03-2026/LSTM_vf_updated_lstm_only_clean_v5_populated-1.ipynb

---
[cell 1 markdown]
# Dow Jones Index Dataset

The **Dow Jones Index** dataset contains weekly stock market data for 30 major companies in the Dow Jones Industrial Average. Each record represents one week of trading activity for a company, including price movements, trading volume, and dividend information.

Although the dataset provides many features, in this study we focus on **predicting the next week’s closing price** using a Long Short-Term Memory (LSTM) network, a special type of Recurrent Neural Network (RNN).

The dataset covers a limited time span: **6 months of 2011 (January–June)** — which results in about 750 weekly records across 30 stocks.

## Source: **UCI Machine Learning Repository**

- [UCI Dataset Page](https://archive.ics.uci.edu/dataset/312/dow+jones+index)

## Target Variable:

- **next_weeks_close** _(numeric)_: Closing price of the stock for the following week.

---

### Features:

| Column Name                        | Description                                                   | Data Type         |
| ---------------------------------- | ------------------------------------------------------------- | ----------------- |
| quarter                            | Fiscal quarter of data (1–4)                                  | _Integer_         |
| stock                              | Stock ticker symbol (e.g., `AA`, `AXP`)                       | _Categorical_     |
| date                               | Week-ending date of the record                                | _Date_            |
| open                               | Stock price at market open (that week)                        | _Numeric (float)_ |
| high                               | Highest stock price during the week                           | _Numeric (float)_ |
| low                                | Lowest stock price during the week                            | _Numeric (float)_ |
| close                              | Stock price at market close (that week)                       | _Numeric (float)_ |
| volume                             | Total trading volume during the week                          | _Numeric (float)_ |
| percent_change_price               | Percent change from opening to closing price                  | _Numeric (float)_ |
| percent_change_volume_over_last_wk | Percent change in volume relative to the previous week        | _Numeric (float)_ |
| previous_weeks_volume              | Trading volume of the previous week                           | _Numeric (float)_ |
| next_weeks_open                    | Stock price at market open (next week)                        | _Numeric (float)_ |
| next_weeks_close                   | Stock price at market close (next week) — **Target**          | **Numeric**       |
| percent_change_next_weeks_price    | Percent change in price from next week’s open to close        | _Numeric (float)_ |
| days_to_next_dividend              | Number of days until the next dividend is paid                | _Integer_         |
| percent_return_next_dividend       | Percent return from the next dividend relative to stock price | _Numeric (float)_ |

[cell 2 markdown]
# **Load into pandas dataframe**

[cell 3 markdown]
# displaying the head

[cell 4 code]
import pandas as pd

# Load file with first row as header
df = pd.read_csv("dow_jones_index.data")

df.head()

[cell 5 code]
df.info()

[cell 6 markdown]
For **predicting next week’s closing price**, we must avoid any feature that already leaks **future information**. In this dataset, the following cannot be used:

- **`next_weeks_open`** → tells us the stock’s opening price in the _future week_.
- **`next_weeks_close`** → is literally the _target variable_, and is redundant since it equals the following row’s `close`, so it cannot be an input feature.
- **`percent_change_next_weeks_price`** → is computed using the next week’s open and close, so it directly encodes future outcome.

[cell 7 markdown]
There are missing values in some columns ( `percent_change_volume_over_last_wk`, `previous_weeks_volume`), but since we are not using these features for our time series modeling, we do not fill the missing values.

[cell 8 markdown]
# **Preprocessing Steps**

- **Clean price columns** : Removed `$` and converted `open`, `high`, `low`, `close`, `next_weeks_open`, `next_weeks_close` to float.
- **Convert date column** : Changed `date` to datetime format for proper time handling.
- **Sort for time series**: Sorted by `stock` and `date` so data is ordered correctly for sequence modeling.

[cell 9 code]
import pandas as pd

def preprocess_dow_jones(path):
    # Load raw CSV
    df = pd.read_csv(path)

    # Remove $ and convert to float for price columns
    money_cols = ["open", "high", "low", "close",
                  "next_weeks_open", "next_weeks_close"]
    for col in money_cols:
        df[col] = df[col].str.replace("$", "", regex=False).astype(float)

    # Convert date column to datetime
    df["date"] = pd.to_datetime(df["date"])

    # Sort values for time-series use (per stock, by date)
    df = df.sort_values(by=["stock", "date"]).reset_index(drop=True)

    return df


df_v1 = preprocess_dow_jones("dow_jones_index.data")
df_v1.head()

[cell 10 code]
df_v1.info()

[cell 11 markdown]
# **Records per Stock**

[cell 12 code]
df_v1["stock"].value_counts()

[cell 13 markdown]
# **Distribution of Closing Prices for Each Stock**

[cell 14 code]
import seaborn as sns
import matplotlib.pyplot as plt

plt.figure(figsize=(12,6))
sns.boxplot(x="stock", y="close", data=df_v1)
plt.xticks(rotation=90)
plt.title("Distribution of Closing Prices by Company")
plt.show()

[cell 15 markdown]
1. Most companies have weekly closing prices that stay within a narrow and stable range i.e. their stock prices do not fluctuate much.
2. A few companies like **IBM** and **CAT** not only have higher average closing prices but also show larger ups and downs, indicating greater volatility.

[cell 16 markdown]
# **Setting fix random seeds across libraries to ensure reproducible results in experiments.**

[cell 17 code]
import random
import numpy as np
import torch

SEED = 42
np.random.seed(SEED)
torch.manual_seed(SEED)
random.seed(SEED)

[cell 18 markdown]
# **Assumption**

[cell 19 markdown]
- **Since each company in the Dow Jones UCI dataset has only about 25 weeks of data, we obtain very few company-specific training examples for LSTM models.**
- **For simplicity and to focus on understanding how time series modeling works, we treat the closing prices in this dataset as if they belonged to a single company’s continuous time series, even though in reality they are grouped by different companies.**

[cell 20 markdown]
# **LSTM Setup for Stock Prediction**

[cell 21 markdown]
# Tuning the Optimal Sequence Length

- The tuning is done with sequence lengths of  4,5, 6,7,8 — meaning we use the closing prices of the past 4,5, 6, 7, or 8 weeks to predict the closing price of the next week, effectively making this a **univariate time series** since only one feature (closing price) is used as input.

[cell 22 markdown]
### **Key Functions & Approach**

- **Creating Sequences**: The core of time-series modeling is transforming the flat list of closing prices into overlapping sequences. The `create_sequences` function creates input windows (`X`) of a specified `seq_len` and corresponding targets (`y`). The choice of `seq_len` is a critical hyperparameter that determines how much past information the model uses.

- **Train–Test Split**: The data is split chronologically (80% for training, 20% for testing) to mimic a real-world scenario where we predict the future based on the past. Random shuffling is avoided as it would destroy the temporal order.

- **Data Scaling**: Stock prices are scaled to a [0, 1] range using **MinMaxScaler**. This helps stabilize the training process and is suitable for price data, which is always non-negative.

- **LSTM Model**: We use an LSTM network, which is well-suited for capturing temporal dependencies. Unlike a simple RNN, an LSTM contains memory cells and gates that allow it to learn and remember information over longer sequences, making it more robust against issues like the vanishing gradient problem.

[cell 23 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import matplotlib.pyplot as plt

# Prepare the data
data = df_v1["close"].values.reshape(-1, 1)
scaler = MinMaxScaler(feature_range=(0, 1))
data_scaled = scaler.fit_transform(data)


torch.manual_seed(SEED)

# Function to create sequences
def create_sequences(data, seq_len=5):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X), np.array(y)

# LSTM Model Definition
class LSTMModel(nn.Module):
    def __init__(self, input_size=1, hidden_size=32, num_layers=1):
        super(LSTMModel, self).__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        # x shape: (batch_size, seq_len, input_size)
        # Get output, final hidden state (hn), and final cell state (cn)
        # We only need the final hidden state
        _, (hn, cn) = self.lstm(x)

        # hn shape: (num_layers, batch_size, hidden_size)
        # We use the hidden state of the last layer for prediction
        out = self.fc(hn[-1])
        return out

# Training and Evaluation Function
def train_and_evaluate(seq_len):
    X, y = create_sequences(data_scaled, seq_len)
    train_size = int(len(X) * 0.8)
    X_train, X_test = X[:train_size], X[train_size:]
    y_train, y_test = y[:train_size], y[train_size:]

    X_train = torch.tensor(X_train, dtype=torch.float32)
    y_train = torch.tensor(y_train, dtype=torch.float32)
    X_test = torch.tensor(X_test, dtype=torch.float32)
    y_test = torch.tensor(y_test, dtype=torch.float32)

    model = LSTMModel()
    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=0.01)

    for epoch in range(50):
        model.train()
        optimizer.zero_grad()
        output = model(X_train)
        loss = criterion(output, y_train)
        loss.backward()
        optimizer.step()

    # Evaluate
    model.eval()
    with torch.no_grad():
        y_pred = model(X_test).numpy()
        y_true = y_test.numpy()

    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    return mse, rmse, mae, r2

# Run experiments for different sequence lengths (NEW RANGE)
seq_lengths = [4,5,6, 7, 8]
results = {"seq_len": [], "mse": [], "rmse": [], "mae": [], "r2": []}

for sl in seq_lengths:
    mse, rmse, mae, r2 = train_and_evaluate(sl)
    results["seq_len"].append(sl)
    results["mse"].append(mse)
    results["rmse"].append(rmse)
    results["mae"].append(mae)
    results["r2"].append(r2)
    print(f"Done: seq_len={sl}")

# Plotting function
def plot_results(results):
    fig, axs = plt.subplots(2, 2, figsize=(12,8))

    axs[0,0].plot(results["seq_len"], results["mse"], marker="o")
    axs[0,0].set_title("MSE vs Seq Length")

    axs[0,1].plot(results["seq_len"], results["rmse"], marker="o")
    axs[0,1].set_title("RMSE vs Seq Length")

    axs[1,0].plot(results["seq_len"], results["mae"], marker="o")
    axs[1,0].set_title("MAE vs Seq Length")

    axs[1,1].plot(results["seq_len"], results["r2"], marker="o")
    axs[1,1].set_title("R² vs Seq Length")

    for ax in axs.flat:
        ax.set_xlabel("Sequence Length")
        ax.grid(True)

    plt.tight_layout()
    plt.show()

plot_results(results)

[cell 24 code]
pd.DataFrame(results)

[cell 25 markdown]
From both the graphs and the tabular results, it can be seen that a sequence length of **7** gives the best performance, achieving the lowest error metrics (MSE, RMSE, MAE) and the highest R² score. This suggests that using the past 6 weeks of data provides the optimal amount of context for the LSTM to predict the next week's price.

[cell 26 markdown]
# **Tuning Hidden Size and Number of Layers in LSTM**

- **Hidden Size**: This hyperparameter controls the number of hidden units (or memory cells) in each LSTM layer. It defines the model's capacity to learn complex patterns. A larger hidden size can capture more intricate dependencies but also increases the risk of overfitting.

- **Number of Layers**: This determines how many LSTM layers are stacked. Stacking layers allows the model to learn hierarchical temporal representations, where each layer processes the output sequence of the previous one. Deeper models can capture more abstract features but are slower to train.

[cell 27 markdown]
After fixing `seq_len` at 7, we performed a grid search over
`hidden_sizes = [32, 64, 128,256]` and `num_layers = [1, 2,3]`,
resulting in a total of 12 combinations of hidden size and number of layers.

[cell 28 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import matplotlib.pyplot as plt

torch.manual_seed(SEED)

# Prepare the data
data = df_v1["close"].values.reshape(-1, 1)
scaler = MinMaxScaler(feature_range=(0, 1))
data_scaled = scaler.fit_transform(data)

def create_sequences(data, seq_len=6):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X), np.array(y)

# LSTM Model Definition
class LSTMModel(nn.Module):
    def __init__(self, input_size=1, hidden_size=32, num_layers=1):
        super(LSTMModel, self).__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (hn, cn) = self.lstm(x)
        # Use the hidden state of the last layer for prediction
        out = self.fc(hn[-1])
        return out

# Training and Evaluation Function
def train_and_evaluate(seq_len, hidden_size, num_layers):
    X, y = create_sequences(data_scaled, seq_len)
    train_size = int(len(X) * 0.8)
    X_train, X_test = X[:train_size], X[train_size:]
    y_train, y_test = y[:train_size], y[train_size:]

    X_train = torch.tensor(X_train, dtype=torch.float32)
    y_train = torch.tensor(y_train, dtype=torch.float32)
    X_test = torch.tensor(X_test, dtype=torch.float32)
    y_test = torch.tensor(y_test, dtype=torch.float32)

    model = LSTMModel(hidden_size=hidden_size, num_layers=num_layers)
    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=0.01)

    for epoch in range(50):
        model.train()
        optimizer.zero_grad()
        output = model(X_train)
        loss = criterion(output, y_train)
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        y_pred = model(X_test).numpy()
        y_true = y_test.numpy()

    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    return mse, rmse, mae, r2

# Grid Search
seq_len = 7 #
hidden_sizes = [32, 64, 128, 256]
num_layers_list = [1, 2, 3]

results = {"hidden_size": [], "num_layers": [], "mse": [], "rmse": [], "mae": [], "r2": []}

for hs in hidden_sizes:
    for nl in num_layers_list:
        mse, rmse, mae, r2 = train_and_evaluate(seq_len, hs, nl)
        results["hidden_size"].append(hs)
        results["num_layers"].append(nl)
        results["mse"].append(mse)
        results["rmse"].append(rmse)
        results["mae"].append(mae)
        results["r2"].append(r2)
        print(f"Done: hidden_size={hs}, num_layers={nl}")

#Plotting
def plot_tuning_results(df_results):
    fig, axs = plt.subplots(2, 2, figsize=(12,8))

    for nl in df_results["num_layers"].unique():
        subset = df_results[df_results["num_layers"]==nl]
        axs[0,0].plot(subset["hidden_size"], subset["mse"], marker="o", label=f"layers={nl}")
        axs[0,1].plot(subset["hidden_size"], subset["rmse"], marker="o", label=f"layers={nl}")
        axs[1,0].plot(subset["hidden_size"], subset["mae"], marker="o", label=f"layers={nl}")
        axs[1,1].plot(subset["hidden_size"], subset["r2"], marker="o", label=f"layers={nl}")

    axs[0,0].set_title("MSE vs Hidden Size")
    axs[0,1].set_title("RMSE vs Hidden Size")
    axs[1,0].set_title("MAE vs Hidden Size")
    axs[1,1].set_title("R² vs Hidden Size")

    for ax in axs.flat:
        ax.set_xlabel("Hidden Size")
        ax.legend()
        ax.grid(True)

    plt.tight_layout()
    plt.show()

df_results = pd.DataFrame(results)
plot_tuning_results(df_results)

[cell 29 code]
df_results

[cell 30 markdown]
From both the graph and the tabular results, a hidden size of **64** with **2 hidden layers** gives the best performance. While a single layer performs well, stacking a second layer with a larger hidden size appears to capture more complex patterns, leading to lower error and a higher R² score. This combination provides the best balance between model capacity and generalization for this dataset.

[cell 31 markdown]
# **Final Model Training and Evaluation**

Based on the hyperparameter tuning, we now train the final LSTM model using the optimal configuration:
- **Sequence Length**: 7
- **Hidden Size**: 64
- **Number of Layers**: 2

We train for 100 epochs to ensure the model converges properly and then evaluate its performance on both the training and test sets.

[cell 32 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler
import matplotlib.pyplot as plt

torch.manual_seed(SEED)

# data prep

data = df_v1["close"].values.reshape(-1, 1)




# Normalize 0-1 scaling
scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(data)

def create_sequences(data, seq_len=6):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X), np.array(y)

# Use optimal sequence length (ADJUST IF YOURS IS DIFFERENT)
SEQ_LEN = 7
X, y = create_sequences(data_scaled, SEQ_LEN)

# Train-test split (chronological)
train_size = int(len(X) * 0.8)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

# Convert to torch tensors
X_train = torch.tensor(X_train, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_test = torch.tensor(y_test, dtype=torch.float32)

# LSTM Model
class LSTMModel(nn.Module):
    def __init__(self, input_size=1, hidden_size=32, num_layers=2):
        super(LSTMModel, self).__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        # The LSTM returns the output sequence and the final hidden and cell states.
        # We only need the final hidden state of the last layer to make a prediction.
        _, (hn, cn) = self.lstm(x)
        # hn is of shape (num_layers, batch_size, hidden_size)
        out = self.fc(hn[-1])
        return out

# Instantiate with optimal hyperparameters (ADJUST IF YOURS ARE DIFFERENT)
model = LSTMModel(hidden_size=128, num_layers=2)

# Train the model
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

EPOCHS = 100
for epoch in range(EPOCHS):
    model.train()
    optimizer.zero_grad()
    output = model(X_train)
    loss = criterion(output, y_train)
    loss.backward()
    optimizer.step()
    if (epoch+1) % 10 == 0:
        print(f"Epoch {epoch+1}/{EPOCHS}, Loss: {loss.item():.6f}")

# Evaluate the model
model.eval()
with torch.no_grad():
    train_pred = model(X_train).numpy()
    test_pred = model(X_test).numpy()

from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

def evaluate_and_print(y_train_true, y_train_pred, y_test_true, y_test_pred):
    def _metrics(y_true, y_pred):
        mse = mean_squared_error(y_true, y_pred)
        rmse = np.sqrt(mse)
        r2 = r2_score(y_true, y_pred)
        return mse, rmse, r2

    train_mse, train_rmse, train_r2 = _metrics(y_train_true, y_train_pred)
    print(f"Train MSE: {train_mse:.4f}, RMSE: {train_rmse:.4f}, R2 score: {train_r2:.4f}")

    test_mse, test_rmse, test_r2 = _metrics(y_test_true, y_test_pred)
    print(f"Test MSE: {test_mse:.4f}, RMSE: {test_rmse:.4f}, R2 score: {test_r2:.4f}")

# Inverse transform to get actual prices
train_pred = scaler.inverse_transform(train_pred)
y_train_real = scaler.inverse_transform(y_train.numpy())
test_pred = scaler.inverse_transform(test_pred)
y_test_real = scaler.inverse_transform(y_test.numpy())

evaluate_and_print(y_train_real, train_pred, y_test_real, test_pred)

[cell 33 markdown]
# **Actual vs Predicted Closing Prices**

[cell 34 code]
plt.figure(figsize=(12,6))
plt.plot(y_test_real, label="Actual Close Price", color="blue")
plt.plot(test_pred, label="Predicted Close Price", color="red", linestyle="--")
plt.title("LSTM: Actual vs. Predicted Closing Prices on Test Set")
plt.legend()
plt.show()

[cell 35 markdown]
**Observations:**

1. The predicted curve follows the actual trend closely, capturing both upward and downward movements with high fidelity.
2. The LSTM model successfully tracks the overall patterns, and even near sharp jumps or drops, the deviations are minimal, demonstrating its strong predictive capability.

[cell 36 markdown]
# **Comparison with MLP**

To benchmark the LSTM's performance, we also trained a simple feedforward neural network (MLP) on the same time-series sequences. The sequences were flattened to be used as input features for the MLP. To highlight the LSTM's ability to handle temporal data, we deliberately use a smaller MLP with a hidden size of 16.

[cell 37 code]
import torch
import torch.nn as nn
import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error

torch.manual_seed(SEED)

scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(df_v1["close"].values.reshape(-1, 1))

def create_sequences(data, seq_len=6):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X), np.array(y)

SEQ_LEN = 7
X, y = create_sequences(data_scaled, SEQ_LEN)

train_size = int(len(X) * 0.8)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

# MLP

class FFN(nn.Module):
    def __init__(self, input_size, hidden_size=16):
        super(FFN, self).__init__()
        self.fc1 = nn.Linear(input_size, hidden_size)
        self.relu = nn.ReLU()
        self.fc2 = nn.Linear(hidden_size, 1)

    def forward(self, x):
        return self.fc2(self.relu(self.fc1(x)))

def train_ffn(X_train, y_train, X_test, y_test, epochs=100):
    input_size = X_train.shape[1]
    model = FFN(input_size=input_size, hidden_size=16)
    criterion = nn.MSELoss()
    optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

    X_train_t = torch.tensor(X_train, dtype=torch.float32)
    y_train_t = torch.tensor(y_train, dtype=torch.float32)
    X_test_t = torch.tensor(X_test, dtype=torch.float32)

    for epoch in range(epochs):
        model.train()
        optimizer.zero_grad()
        output = model(X_train_t)
        loss = criterion(output, y_train_t)
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        train_pred = model(X_train_t).numpy()
        test_pred = model(X_test_t).numpy()
    return train_pred, test_pred

X_train_flat = X_train.reshape(X_train.shape[0], -1)
X_test_flat = X_test.reshape(X_test.shape[0], -1)
train_pred_ffn, test_pred_ffn = train_ffn(X_train_flat, y_train, X_test_flat, y_test)

# inverse transform for MLP
train_pred_ffn = scaler.inverse_transform(train_pred_ffn)
test_pred_ffn = scaler.inverse_transform(test_pred_ffn)
y_train_real_ffn = scaler.inverse_transform(y_train)
y_test_real_ffn = scaler.inverse_transform(y_test)

# RNN

class RNNModel(nn.Module):
    def __init__(self, input_size=1, hidden_size=32, num_layers=1):
        super(RNNModel, self).__init__()
        self.rnn = nn.RNN(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        out, h = self.rnn(x)
        out = self.fc(h[-1])
        return out

def train_rnn(X_train, y_train, X_test, y_test, epochs=100):
    model = RNNModel(input_size=1, hidden_size=32, num_layers=1)
    criterion = nn.MSELoss()
    optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

    X_train_t = torch.tensor(X_train, dtype=torch.float32)
    y_train_t = torch.tensor(y_train, dtype=torch.float32)
    X_test_t = torch.tensor(X_test, dtype=torch.float32)

    for epoch in range(epochs):
        model.train()
        optimizer.zero_grad()
        output = model(X_train_t)
        loss = criterion(output, y_train_t)
        loss.backward()
        optimizer.step()

    model.eval()
    with torch.no_grad():
        train_pred = model(X_train_t).numpy()
        test_pred = model(X_test_t).numpy()
    return train_pred, test_pred

train_pred_rnn, test_pred_rnn = train_rnn(X_train, y_train, X_test, y_test)

# inverse transform for RNN
train_pred_rnn = scaler.inverse_transform(train_pred_rnn)
test_pred_rnn = scaler.inverse_transform(test_pred_rnn)
y_train_real_rnn = scaler.inverse_transform(y_train)
y_test_real_rnn = scaler.inverse_transform(y_test)

# Evaluation

def evaluate(y_train, train_pred, y_test, test_pred, name):
    def _metrics(y_true, y_pred):
        mse = mean_squared_error(y_true, y_pred)
        rmse = np.sqrt(mse)
        mae = mean_absolute_error(y_true, y_pred)
        r2 = r2_score(y_true, y_pred)
        return mse, rmse, mae, r2

    train_mse, train_rmse, train_mae, train_r2 = _metrics(y_train, train_pred)
    test_mse, test_rmse, test_mae, test_r2 = _metrics(y_test, test_pred)

    return {
        "Model": name,
        "Train MSE": train_mse, "Train RMSE": train_rmse, "Train MAE": train_mae, "Train R2": train_r2,
        "Test MSE": test_mse, "Test RMSE": test_rmse, "Test MAE": test_mae, "Test R2": test_r2
    }

results_list = []
results_list.append(evaluate(y_train_real_ffn, train_pred_ffn, y_test_real_ffn, test_pred_ffn, "MLP"))
results_list.append(evaluate(y_train_real_rnn, train_pred_rnn, y_test_real_rnn, test_pred_rnn, "RNN"))
results_list.append(evaluate(y_train_real, train_pred, y_test_real, test_pred, "LSTM"))

df_results = pd.DataFrame(results_list)
df_results[['Model','Test MSE','Test RMSE','Test MAE','Test R2']]

[cell 38 markdown]
The results show that the performance of LSTM and RNN is comparable and they outperforms the MLP. This highlights the LSTM's/RNN strength in modeling sequential data, as it can capture temporal dependencies that a simple MLP with flattened features misses.

[cell 39 markdown]
# **Multi-Step Forecasting: Predicting the Next 4 Weeks**

Here, the LSTM is adapted to perform **multi-step forecasting**. Instead of predicting just the next week’s price, the model will output predictions for the next 4 weeks simultaneously. This is a sequence-to-sequence (Seq2Seq) task where the input is a sequence of past prices and the output is a sequence of future prices.

The model uses an encoder-decoder architecture. The encoder (an LSTM) processes the input sequence and summarizes it into a context vector (the final hidden and cell states). The decoder (another LSTM) then uses this context to generate the output sequence one step at a time, feeding its own previous prediction back as input for the next step (autoregression).

[cell 40 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error

# Prepare data
data = df_v1["close"].values.reshape(-1, 1)
scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(data)

# Create sequences for multi-step prediction
def create_sequences_multi_step(data, seq_len=6, pred_horizon=4):
    X, y = [], []
    for i in range(len(data) - seq_len - pred_horizon + 1):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len : i+seq_len+pred_horizon].flatten())
    return np.array(X), np.array(y)

SEQ_LEN = 7
PRED_HORIZON = 4
X, y = create_sequences_multi_step(data_scaled, SEQ_LEN, PRED_HORIZON)

train_size = int(len(X) * 0.8)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

X_train = torch.tensor(X_train, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_test = torch.tensor(y_test, dtype=torch.float32)

# Seq2Seq LSTM model
class Seq2SeqLSTM(nn.Module):
    def __init__(self, input_size=1, hidden_size=32, num_layers=2, pred_horizon=4):
        super(Seq2SeqLSTM, self).__init__()
        self.pred_horizon = pred_horizon
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x, y=None, train=False):
        # Encoder
        _, (h, c) = self.lstm(x)

        # Start decoding with the last observed value
        dec_input = x[:, -1:, :]   # (batch, 1, input_size)
        preds = []

        for t in range(self.pred_horizon):
            out, (h, c) = self.lstm(dec_input, (h, c))
            y_pred = self.fc(out[:, -1, :])  # (batch, 1)
            preds.append(y_pred)

            if train and y is not None:
                # ground truth as next input during training
                dec_input = y[:, t].unsqueeze(1).unsqueeze(2)  # (batch, 1, 1)
            else:
                # own prediction during inference
                dec_input = y_pred.unsqueeze(1)

        return torch.cat(preds, dim=1)  # (batch, pred_horizon)


model = Seq2SeqLSTM(hidden_size=64, num_layers=2, pred_horizon=PRED_HORIZON)

# Training
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

EPOCHS = 100
for epoch in range(EPOCHS):
    model.train()
    optimizer.zero_grad()
    output = model(X_train, y_train, train=True)
    loss = criterion(output, y_train)
    loss.backward()
    optimizer.step()
    if (epoch+1) % 10 == 0:
        print(f"Epoch {epoch+1}/{EPOCHS}, Loss: {loss.item():.6f}")

# Evaluation
model.eval()
with torch.no_grad():
    test_pred = model(X_test, train=False).numpy()

# Inverse transform
test_pred = scaler.inverse_transform(test_pred)
y_test_real = scaler.inverse_transform(y_test.numpy())


# Metrics
def evaluate_multi_step(y_true, y_pred):
    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    mae = mean_absolute_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    print(f"Test MSE: {mse:.4f}, RMSE: {rmse:.4f}, MAE: {mae:.4f}, R2: {r2:.4f}")

evaluate_multi_step(y_test_real, test_pred)

[cell 41 code]
import matplotlib.pyplot as plt

fig, axs = plt.subplots(2, 2, figsize=(14,8))

# Week +1
axs[0,0].plot(y_test_real[:,0], label="Actual Week+1")
axs[0,0].plot(test_pred[:,0], label="Predicted Week+1", linestyle="--")
axs[0,0].legend()
axs[0,0].set_title("Week +1")

# Week +2
axs[0,1].plot(y_test_real[:,1], label="Actual Week+2")
axs[0,1].plot(test_pred[:,1], label="Predicted Week+2", linestyle="--")
axs[0,1].legend()
axs[0,1].set_title("Week +2")

# Week +3
axs[1,0].plot(y_test_real[:,2], label="Actual Week+3")
axs[1,0].plot(test_pred[:,2], label="Predicted Week+3", linestyle="--")
axs[1,0].legend()
axs[1,0].set_title("Week +3")

# Week +4
axs[1,1].plot(y_test_real[:,3], label="Actual Week+4")
axs[1,1].plot(test_pred[:,3], label="Predicted Week+4", linestyle="--")
axs[1,1].legend()
axs[1,1].set_title("Week +4")

plt.suptitle("Multi-step LSTM Prediction (Next 4 Weeks Closing Price)", fontsize=14)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()

[cell 42 markdown]
The performance clearly dropped in the multi-step setting compared to single-step prediction. This is expected, as predicting further into the future is inherently more difficult. Errors accumulate in the autoregressive decoding process: a small error in the Week +1 prediction is fed as input to predict Week +2, potentially causing a larger error, and so on.

Despite this, the model captures the overall trends well for all four weeks. As the prediction horizon increases from Week +1 to Week +4, the deviation from the actual values becomes more pronounced, but the general direction is often correct.

[cell 43 markdown]
# **LSTM with Multiple Features for Prediction**

In this final experiment, we enhance the LSTM model by using multiple input features instead of just the closing price. We use a multivariate time series including `open`, `high`, `low`, `close`, and `volume` (OHLCV) from the past 6 weeks to predict the next week’s closing price. This allows the model to learn from richer, more contextual information about the stock's weekly behavior.

[cell 44 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error

torch.manual_seed(SEED)

feature_cols = ["open", "high", "low", "close", "volume"]
data_multi = df_v1[feature_cols].values
target = df_v1["close"].values.reshape(-1, 1)

# Normalize features and target separately
scaler_X = MinMaxScaler()
scaler_y = MinMaxScaler()
data_scaled = scaler_X.fit_transform(data_multi)
target_scaled = scaler_y.fit_transform(target)

SEQ_LEN = 7 # Use the optimal sequence length from tuning

def create_sequences_multi_feature(X, y, seq_len=6):
    Xs, ys = [], []
    for i in range(len(X) - seq_len):
        Xs.append(X[i:i+seq_len])
        ys.append(y[i+seq_len])
    return np.array(Xs), np.array(ys)

X, y = create_sequences_multi_feature(data_scaled, target_scaled, SEQ_LEN)
print("X shape:", X.shape)
print("y shape:", y.shape)

# Train-test split
train_size = int(len(X) * 0.8)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

X_train = torch.tensor(X_train, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_test = torch.tensor(y_test, dtype=torch.float32)

# LSTM model (multi-feature)
class LSTMModelMulti(nn.Module):
    def __init__(self, input_size, hidden_size=32, num_layers=1):
        super(LSTMModelMulti, self).__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (hn, cn) = self.lstm(x)
        # Use the final hidden state from the last layer for prediction
        out = self.fc(hn[-1])
        return out

num_features = X.shape[2]
# Use optimal hyperparameters
model = LSTMModelMulti(input_size=num_features, hidden_size=128, num_layers=2)

# Training
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

EPOCHS = 100
for epoch in range(EPOCHS):
    model.train()
    optimizer.zero_grad()
    output = model(X_train)
    loss = criterion(output, y_train)
    loss.backward()
    optimizer.step()
    if (epoch+1) % 10 == 0:
        print(f"Epoch {epoch+1}/{EPOCHS}, Loss: {loss.item():.6f}")

# Evaluation
model.eval()
with torch.no_grad():
    train_pred = model(X_train).numpy()
    test_pred = model(X_test).numpy()

# Inverse transform
train_pred = scaler_y.inverse_transform(train_pred)
y_train_real = scaler_y.inverse_transform(y_train.numpy())
test_pred = scaler_y.inverse_transform(test_pred)
y_test_real = scaler_y.inverse_transform(y_test.numpy())

def evaluate_and_print(y_train_true, y_train_pred, y_test_true, y_test_pred):
    def _metrics(y_true, y_pred):
        mse = mean_squared_error(y_true, y_pred)
        rmse = np.sqrt(mse)
        r2 = r2_score(y_true, y_pred)
        mae=mean_absolute_error(y_true,y_pred)
        return mse, rmse, r2,mae

    train_mse, train_rmse, train_r2,train_mae = _metrics(y_train_true, y_train_pred)
    print(f"Train MSE: {train_mse:.4f}, RMSE: {train_rmse:.4f}, MAE: {train_mae:.4f}, R2 score: {train_r2:.4f}")

    test_mse, test_rmse, test_r2,test_mae = _metrics(y_test_true, y_test_pred)
    print(f"Test MSE: {test_mse:.4f}, RMSE: {test_rmse:.4f}, MAE: {test_mae:.4f}, R2 score: {test_r2:.4f}")


evaluate_and_print(y_train_real, train_pred, y_test_real, test_pred)

[cell 45 markdown]
Using multiple input features led to a slight dip in performance, suggesting that additional variables may be redundant and the closing price alone already captures most of the predictive signal.

[cell 46 markdown]
# **Implementing LSTM from Scratch**

An LSTM cell controls information flow using four gates:

- **Forget Gate ($f_t$)** → controls how much of the previous cell state ($c_{t-1}$) to keep.  
- **Input Gate ($i_t$)** → decides how much new candidate information should enter the cell.  
- **Candidate State ($g_t$)** → new content that can be added into the cell state.  
- **Output Gate ($o_t$)** → decides how much of the cell state influences the hidden state ($h_t$).  

### Updates
$$
c_t = f_t \odot c_{t-1} + i_t \odot g_t
$$

$$
h_t = o_t \odot \tanh(c_t)
$$

[cell 47 code]
import torch
import torch.nn as nn
import math

class SimpleLSTM(nn.Module):
    def __init__(self, input_size, hidden_size):
        super(SimpleLSTM, self).__init__()
        self.input_size = input_size
        self.hidden_size = hidden_size

        # One weight matrix for input -> gates
        self.W_x = nn.Parameter(torch.Tensor(input_size, 4 * hidden_size))
        # One weight matrix for hidden -> gates
        self.W_h = nn.Parameter(torch.Tensor(hidden_size, 4 * hidden_size))
        # One bias vector
        self.b = nn.Parameter(torch.Tensor(4 * hidden_size))

        self.reset_parameters()

    def reset_parameters(self):
        # Initialize weights uniformly (like PyTorch does)
        stdv = 1.0 / math.sqrt(self.hidden_size)
        for param in self.parameters():
            nn.init.uniform_(param, -stdv, stdv)

    def forward(self, x, state=None):
        """
        x: (batch, seq_len, input_size)
        state: (h0, c0), each (batch, hidden_size)
        """
        batch, seq_len, _ = x.size()

        if state is None:
            h_t = torch.zeros(batch, self.hidden_size, device=x.device)
            c_t = torch.zeros(batch, self.hidden_size, device=x.device)
        else:
            h_t, c_t = state

        outputs = []

        for t in range(seq_len):
            x_t = x[:, t, :]  # (batch, input_size)

            gates = x_t @ self.W_x + h_t @ self.W_h + self.b
            i, f, g, o = gates.chunk(4, dim=1)  # split into 4 parts

            i = torch.sigmoid(i)      # input gate
            f = torch.sigmoid(f)      # forget gate
            g = torch.tanh(g)         # candidate cell
            o = torch.sigmoid(o)      # output gate

            c_t = f * c_t + i * g     # new cell state
            h_t = o * torch.tanh(c_t) # new hidden state

            outputs.append(h_t.unsqueeze(1))

        outputs = torch.cat(outputs, dim=1)  # (batch, seq_len, hidden_size)
        return outputs, (h_t, c_t)


#  Example
lstm = SimpleLSTM(input_size=10, hidden_size=20)

x = torch.randn(3, 5, 10)   # (batch=3, seq_len=5, input_size=10)
h0 = torch.randn(3, 20)
c0 = torch.randn(3, 20)

out, (hn, cn) = lstm(x, (h0, c0))

print(out.shape)  # (3, 5, 20)
print(hn.shape)   # (3, 20)
print(cn.shape)   # (3, 20)

[cell 48 markdown]
# **Stacked LSTM**

- **Manual stacking:** just pass the output of one LSTM as the input to the next.

[cell 49 code]
inp = torch.randn(3, 5, 10)

lstm1 = SimpleLSTM(10, 20)
lstm2 = SimpleLSTM(20, 20)

out1, (h1, c1) = lstm1(inp)
out2, (h2, c2) = lstm2(out1)

print(out1.shape)
print(out2.shape)

[cell 50 markdown]
- **Define a `StackedLSTM` class**. It creates multiple `SimpleLSTM` layers in a loop and stores them in a `ModuleList`. The forward pass also loops through these layers. This makes it easy to stack many layers by just setting `num_layers`.

[cell 51 code]
class StackedLSTM(nn.Module):
    def __init__(self, input_size, hidden_size, num_layers):
        super(StackedLSTM, self).__init__()
        layers = []
        for i in range(num_layers):
            in_size = input_size if i == 0 else hidden_size
            layers.append(SimpleLSTM(in_size, hidden_size))
        self.layers = nn.ModuleList(layers)

    def forward(self, x):
        """
        x: (batch, seq_len, input_size)
        """
        out = x
        for layer in self.layers:
            out, _ = layer(out)   # each SimpleLSTM keeps batch-first
        return out


# Example
inp = torch.randn(3, 5, 10)   # (batch=3, seq_len=5, input_size=10)
stacked = StackedLSTM(10, 20, num_layers=2)
out = stacked(inp)
print(out.shape)   # (3, 5, 20)

[cell 52 markdown]
# **LSTM Forecaster with manually implemented LSTM (One-Step Prediction)**
LSTM with a linear head to do prediction

[cell 53 code]
import torch
import torch.nn as nn
import math

class LSTMForecaster(nn.Module):
    def __init__(self, input_size, hidden_size, output_size):
        super(LSTMForecaster, self).__init__()
        self.lstm = SimpleLSTM(input_size, hidden_size)
        self.fc = nn.Linear(hidden_size, output_size)  # prediction head

    def forward(self, x, state=None):
        """
        x: (batch, seq_len, input_size)
        """
        outputs, (hn, cn) = self.lstm(x, state)
        # Use last hidden state for one-step prediction
        pred = self.fc(hn)   # hn shape: (batch, hidden_size)
        return pred


# ---------------- Example ----------------
torch.manual_seed(0)

model = LSTMForecaster(input_size=10, hidden_size=20, output_size=1)

# Dummy input: batch=3, seq_len=5, input_size=10
x = torch.randn(3, 5, 10)

# Dummy target: one-step prediction -> shape (batch, output_size)
target = torch.randn(3, 1)

# Forward pass
pred = model(x)
print("Prediction shape:", pred.shape)   # (3, 1)

# Loss
criterion = nn.MSELoss()
loss = criterion(pred, target)
print("Loss:", loss.item())

# Backward pass
loss.backward()

# Print gradients
print("\nGradients:")
print("W_x.grad norm:", model.lstm.W_x.grad.norm().item())
print("W_h.grad norm:", model.lstm.W_h.grad.norm().item())
print("b.grad norm:", model.lstm.b.grad.norm().item())
print("fc.weight.grad norm:", model.fc.weight.grad.norm().item())
print("fc.bias.grad norm:", model.fc.bias.grad.norm().item())

[cell 54 markdown]
# **Using manually implemented LSTM to do one step closing price prediction**

[cell 55 code]
import pandas as pd
import numpy as np
import torch
import torch.nn as nn
from sklearn.preprocessing import MinMaxScaler
import matplotlib.pyplot as plt

data = df_v1["close"].values.reshape(-1, 1)

torch.manual_seed(SEED)

# Normalize 0-1 scaling
scaler = MinMaxScaler()
data_scaled = scaler.fit_transform(data)

def create_sequences(data, seq_len=6):
    X, y = [], []
    for i in range(len(data) - seq_len):
        X.append(data[i:i+seq_len])
        y.append(data[i+seq_len])
    return np.array(X), np.array(y)

# Use optimal sequence length (ADJUST IF YOURS IS DIFFERENT)
SEQ_LEN = 7
X, y = create_sequences(data_scaled, SEQ_LEN)

print(X.shape,y.shape)

# Train-test split (chronological)
train_size = int(len(X) * 0.8)
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

# Convert to torch tensors
X_train = torch.tensor(X_train, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_test = torch.tensor(y_test, dtype=torch.float32)


# Instantiate with optimal hyperparameters (ADJUST IF YOURS ARE DIFFERENT)
model = LSTMForecaster(input_size=1, hidden_size=128, output_size=1)

# Train the model
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.01)

EPOCHS = 100
for epoch in range(EPOCHS):
    model.train()
    optimizer.zero_grad()
    output = model(X_train)
    #print(X_train.shape)
    #print(output.shape,y_train.shape)
    loss = criterion(output, y_train)
    loss.backward()
    optimizer.step()
    if (epoch+1) % 10 == 0:
        print(f"Epoch {epoch+1}/{EPOCHS}, Loss: {loss.item():.6f}")

# Evaluate the model
model.eval()
with torch.no_grad():
    train_pred = model(X_train).numpy()
    test_pred = model(X_test).numpy()

from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

def evaluate_and_print(y_train_true, y_train_pred, y_test_true, y_test_pred):
    def _metrics(y_true, y_pred):
        mse = mean_squared_error(y_true, y_pred)
        rmse = np.sqrt(mse)
        r2 = r2_score(y_true, y_pred)
        return mse, rmse, r2

    train_mse, train_rmse, train_r2 = _metrics(y_train_true, y_train_pred)
    print(f"Train MSE: {train_mse:.4f}, RMSE: {train_rmse:.4f}, R2 score: {train_r2:.4f}")

    test_mse, test_rmse, test_r2 = _metrics(y_test_true, y_test_pred)
    print(f"Test MSE: {test_mse:.4f}, RMSE: {test_rmse:.4f}, R2 score: {test_r2:.4f}")


evaluate_and_print(y_train, train_pred, y_test, test_pred)

[cell 56 markdown]
# **LSTM for Classification: EEG Eye State**

Here we use the EEG Eye State dataset to perform binary classification of eye state (open vs. closed) from EEG signals. This section compares the performance of an LSTM against a standard MLP to see if the temporal nature of the EEG signal provides an advantage.

[cell 57 code]
!pip install -q liac-arff

[cell 58 code]
# Download the EEG Eye State dataset from the UCI repository
!wget https://archive.ics.uci.edu/ml/machine-learning-databases/00264/EEG%20Eye%20State.arff

[cell 59 code]
from scipy.io import arff
import pandas as pd

# Load ARFF
data, meta = arff.loadarff("EEG Eye State.arff")
df_eeg = pd.DataFrame(data)

# Decode byte strings and convert target to integer
df_eeg['eyeDetection'] = df_eeg['eyeDetection'].str.decode('utf-8').astype(int)
df_eeg.head()

[cell 60 markdown]
# **Training a MLP for Comparison**

[cell 61 code]
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# Prepare data
X = df_eeg.drop('eyeDetection', axis=1).astype(float).values
y = df_eeg['eyeDetection'].astype(int).values

# Split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# Scale
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Torch tensors
X_train_t = torch.tensor(X_train_scaled, dtype=torch.float32)
y_train_t = torch.tensor(y_train, dtype=torch.long)
X_test_t  = torch.tensor(X_test_scaled, dtype=torch.float32)
y_test_t  = torch.tensor(y_test, dtype=torch.long)

train_ds = TensorDataset(X_train_t, y_train_t)
train_loader = DataLoader(train_ds, batch_size=64, shuffle=True)

# MLP Model
class FFN(nn.Module):
    def __init__(self, input_size, hidden_size=64):
        super(FFN, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_size, hidden_size), nn.ReLU(),
            nn.Linear(hidden_size, hidden_size), nn.ReLU(),
            nn.Linear(hidden_size, 2)
        )
    def forward(self, x): return self.net(x)

model = FFN(input_size=X_train_scaled.shape[1])
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# Training loop
for epoch in range(10):
    for xb, yb in train_loader:
        preds = model(xb)
        loss = criterion(preds, yb)
        loss.backward()
        optimizer.step()
        optimizer.zero_grad()
    print(f"Epoch {epoch+1}, Loss: {loss.item():.4f}")

# Evaluation
with torch.no_grad():
    y_pred = model(X_test_t).argmax(dim=1).numpy()
    y_proba = torch.softmax(model(X_test_t), dim=1)[:,1].numpy()

ffn_results = {
    "Model": "MLP",
    "Accuracy": accuracy_score(y_test, y_pred),
    "Precision": precision_score(y_test, y_pred),
    "Recall": recall_score(y_test, y_pred),
    "F1": f1_score(y_test, y_pred),
    "ROC-AUC": roc_auc_score(y_test, y_proba)
}
print(ffn_results)

[cell 62 markdown]
# **Using LSTM for Predicting Eye State**

Here sequences are artificially created from the continuous EEG stream, so that after every 50 timesteps the model predicts the eye state at the next timestep. This framing lets the LSTM learn temporal dependencies instead of treating each reading independently.

[cell 63 code]
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score

# Data already split and scaled from previous cell
# Create sequences
SEQ_LEN = 50

def create_sequences_clf(X, y, seq_len=50):
    Xs, ys = [], []
    for i in range(len(X) - seq_len):
        Xs.append(X[i:i+seq_len])
        ys.append(y[i+seq_len])
    return np.array(Xs), np.array(ys)

X_train_seq, y_train_seq = create_sequences_clf(X_train_scaled, y_train, SEQ_LEN)
X_test_seq, y_test_seq   = create_sequences_clf(X_test_scaled, y_test, SEQ_LEN)

print("LSTM Train shape:", X_train_seq.shape, y_train_seq.shape)
print("LSTM Test shape :", X_test_seq.shape, y_test_seq.shape)

# Torch tensors
X_train_t = torch.tensor(X_train_seq, dtype=torch.float32)
y_train_t = torch.tensor(y_train_seq, dtype=torch.float32)
X_test_t  = torch.tensor(X_test_seq, dtype=torch.float32)
y_test_t  = torch.tensor(y_test_seq, dtype=torch.float32)

train_ds = TensorDataset(X_train_t, y_train_t)
train_loader = DataLoader(train_ds, batch_size=64, shuffle=True)

# LSTM for Binary Classification
class LSTMBinary(nn.Module):
    def __init__(self, input_size, hidden_size=64, num_layers=2):
        super(LSTMBinary, self).__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (h, _) = self.lstm(x)
        out = self.fc(h[-1]) # Use last layer's hidden state
        return out.squeeze(1)

model = LSTMBinary(input_size=X_train_seq.shape[2])
criterion = nn.BCEWithLogitsLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# Training loop
for epoch in range(10):
    for xb, yb in train_loader:
        preds = model(xb)
        loss = criterion(preds, yb)
        loss.backward()
        optimizer.step()
        optimizer.zero_grad()
    print(f"Epoch {epoch+1}, Loss: {loss.item():.4f}")

# Evaluation
with torch.no_grad():
    y_proba = torch.sigmoid(model(X_test_t)).numpy()
    y_pred  = (y_proba >= 0.5).astype(int)

lstm_results = {
    "Model": "LSTM",
    "Accuracy": accuracy_score(y_test_seq, y_pred),
    "Precision": precision_score(y_test_seq, y_pred),
    "Recall": recall_score(y_test_seq, y_pred),
    "F1": f1_score(y_test_seq, y_pred),
    "ROC-AUC": roc_auc_score(y_test_seq, y_proba)
}
print(lstm_results)

[cell 64 code]
import pandas as pd


results_df = pd.DataFrame([ffn_results, lstm_results])

results_df

[cell 65 markdown]
**MLP gave better results, while the LSTM performed poorly. This shows that forcing sequential framing does not help for this dataset and a plain MLP works better.**

[cell 67 markdown]
## Additional LSTM Demo: Stock-aware LSTM with ticker embeddings (PyTorch)

This section trains a single LSTM across all Dow Jones stocks while giving the model a learned embedding for the ticker symbol. It helps demonstrate how LSTMs can share patterns across related sequences while still modeling per-series differences.

[cell 68 code]
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
import matplotlib.pyplot as plt

# Ensure the preprocessed Dow Jones dataframe is available
if "df_v1" not in globals():
    df_v1 = preprocess_dow_jones("dow_jones_index.data")

torch.manual_seed(42)
np.random.seed(42)

FEATURE_COLS = ["open", "high", "low", "close", "volume"]
TARGET_COL = "close"
SEQ_LEN = 8  # different from earlier sections on purpose

# Build per-stock sequences (chronological split per stock)
stocks = sorted(df_v1["stock"].unique())  # ticker symbols (company identifiers) in the dataset
stock_to_id = {s: i for i, s in enumerate(stocks)}  # map ticker -> integer id for embedding lookup
id_to_stock = {i: s for s, i in stock_to_id.items()}  # inverse map: id -> ticker (for readable outputs)
# The embedding is a trainable vector per ticker id that conditions the LSTM on the stock identity.



def build_sequences_for_stock(g: pd.DataFrame, seq_len: int):
    g = g.sort_values("date").reset_index(drop=True)
    X_raw = g[FEATURE_COLS].astype(float).values
    y_raw = g[TARGET_COL].astype(float).values.reshape(-1, 1)
    dates = g["date"].values

    X, y, d = [], [], []
    for i in range(len(g) - seq_len):
        X.append(X_raw[i:i + seq_len])
        y.append(y_raw[i + seq_len])
        d.append(dates[i + seq_len])
    return np.array(X), np.array(y), np.array(d)

X_train_list, y_train_list, sid_train_list, date_train_list = [], [], [], []
X_test_list, y_test_list, sid_test_list, date_test_list = [], [], [], []

for s, g in df_v1.groupby("stock"):
    X_s, y_s, d_s = build_sequences_for_stock(g, SEQ_LEN)
    if len(X_s) < 5:
        continue

    split = int(len(X_s) * 0.8)
    sid = stock_to_id[s]

    X_train_list.append(X_s[:split])
    y_train_list.append(y_s[:split])
    sid_train_list.append(np.full((split,), sid, dtype=np.int64))
    date_train_list.append(d_s[:split])

    X_test_list.append(X_s[split:])
    y_test_list.append(y_s[split:])
    sid_test_list.append(np.full((len(X_s) - split,), sid, dtype=np.int64))
    date_test_list.append(d_s[split:])

X_train_raw = np.concatenate(X_train_list, axis=0)
y_train_raw = np.concatenate(y_train_list, axis=0)
sid_train = np.concatenate(sid_train_list, axis=0)
date_train = np.concatenate(date_train_list, axis=0)

X_test_raw = np.concatenate(X_test_list, axis=0)
y_test_raw = np.concatenate(y_test_list, axis=0)
sid_test = np.concatenate(sid_test_list, axis=0)
date_test = np.concatenate(date_test_list, axis=0)

# Fit scalers on training data only (avoid test leakage)
scaler_X = MinMaxScaler()
scaler_y = MinMaxScaler()

X_train_2d = X_train_raw.reshape(-1, len(FEATURE_COLS))
X_test_2d = X_test_raw.reshape(-1, len(FEATURE_COLS))

X_train_scaled = scaler_X.fit_transform(X_train_2d).reshape(X_train_raw.shape)
X_test_scaled = scaler_X.transform(X_test_2d).reshape(X_test_raw.shape)

y_train_scaled = scaler_y.fit_transform(y_train_raw)
y_test_scaled = scaler_y.transform(y_test_raw)

# Torch tensors
X_train_t = torch.tensor(X_train_scaled, dtype=torch.float32)
y_train_t = torch.tensor(y_train_scaled, dtype=torch.float32)
sid_train_t = torch.tensor(sid_train, dtype=torch.long)

X_test_t = torch.tensor(X_test_scaled, dtype=torch.float32)
y_test_t = torch.tensor(y_test_scaled, dtype=torch.float32)
sid_test_t = torch.tensor(sid_test, dtype=torch.long)

train_loader = DataLoader(TensorDataset(X_train_t, sid_train_t, y_train_t), batch_size=32, shuffle=True)
test_loader = DataLoader(TensorDataset(X_test_t, sid_test_t, y_test_t), batch_size=64, shuffle=False)

print(f"Train samples: {len(X_train_t):,} | Test samples: {len(X_test_t):,} | Stocks: {len(stocks)}")

[cell 69 code]
class StockAwareLSTM(nn.Module):
    def __init__(self, num_features: int, num_stocks: int, emb_dim: int = 8, hidden_size: int = 64, num_layers: int = 2, dropout: float = 0.2):
        super().__init__()
        self.emb = nn.Embedding(num_stocks, emb_dim)
        self.lstm = nn.LSTM(
            input_size=num_features + emb_dim,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=dropout if num_layers > 1 else 0.0
        )
        self.head = nn.Sequential(
            nn.Linear(hidden_size, 32),
            nn.ReLU(),
            nn.Linear(32, 1)
        )

    def forward(self, x: torch.Tensor, stock_id: torch.Tensor) -> torch.Tensor:
        e = self.emb(stock_id)  # (batch, emb_dim)
        e = e.unsqueeze(1).expand(-1, x.size(1), -1)  # repeat for each timestep
        x_in = torch.cat([x, e], dim=-1)
        out, _ = self.lstm(x_in)
        last = out[:, -1, :]
        return self.head(last)

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = StockAwareLSTM(num_features=len(FEATURE_COLS), num_stocks=len(stocks)).to(device)

criterion = nn.MSELoss()
optimizer = optim.Adam(model.parameters(), lr=1e-3)

best_val = float("inf")
best_state = None
patience = 20
stale = 0
EPOCHS = 200

for epoch in range(1, EPOCHS + 1):
    model.train()
    train_losses = []
    for xb, sidb, yb in train_loader:
        xb, sidb, yb = xb.to(device), sidb.to(device), yb.to(device)
        optimizer.zero_grad()
        pred = model(xb, sidb)
        loss = criterion(pred, yb)
        loss.backward()
        optimizer.step()
        train_losses.append(loss.item())

    model.eval()
    val_losses = []
    with torch.no_grad():
        for xb, sidb, yb in test_loader:
            xb, sidb, yb = xb.to(device), sidb.to(device), yb.to(device)
            pred = model(xb, sidb)
            val_losses.append(criterion(pred, yb).item())
    val_loss = float(np.mean(val_losses))

    if val_loss < best_val - 1e-6:
        best_val = val_loss
        best_state = {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}
        stale = 0
    else:
        stale += 1

    if epoch % 20 == 0 or epoch == 1:
        print(f"Epoch {epoch:3d} | train_loss={np.mean(train_losses):.5f} | val_loss={val_loss:.5f}")

    if stale >= patience:
        print(f"Early stopping at epoch {epoch} (best val_loss={best_val:.5f})")
        break

if best_state is not None:
    model.load_state_dict(best_state)
model.eval();

[cell 70 code]
# Evaluate on the held-out time slices (per stock)
all_pred = []
all_true = []
all_sid = []

with torch.no_grad():
    for xb, sidb, yb in test_loader:
        xb, sidb = xb.to(device), sidb.to(device)
        pred = model(xb, sidb).cpu().numpy()
        all_pred.append(pred)
        all_true.append(yb.numpy())
        all_sid.append(sidb.cpu().numpy())

pred_scaled = np.vstack(all_pred)
true_scaled = np.vstack(all_true)
sid_all = np.concatenate(all_sid)

pred_real = scaler_y.inverse_transform(pred_scaled).ravel()
true_real = scaler_y.inverse_transform(true_scaled).ravel()

rmse = float(np.sqrt(mean_squared_error(true_real, pred_real)))
mae = mean_absolute_error(true_real, pred_real)
r2 = r2_score(true_real, pred_real)

print(f"Stock-aware LSTM test metrics | RMSE: {rmse:.4f} | MAE: {mae:.4f} | R2: {r2:.4f}")

pred_df = pd.DataFrame({
    "stock": [id_to_stock[int(i)] for i in sid_all],
    "date": date_test[:len(sid_all)],
    "actual_close": true_real,
    "pred_close": pred_real
}).sort_values(["stock", "date"]).reset_index(drop=True)

pred_df.head()

[cell 71 code]
# Plot a single ticker to visualize performance
ticker = pred_df["stock"].value_counts().index[0]
plot_df = pred_df[pred_df["stock"] == ticker].sort_values("date")

plt.figure(figsize=(10, 4))
plt.plot(plot_df["date"], plot_df["actual_close"], label="Actual")
plt.plot(plot_df["date"], plot_df["pred_close"], label="Predicted", linestyle="--")
plt.title(f"Stock-aware LSTM: next-week close prediction (ticker={ticker})")
plt.xlabel("Date")
plt.ylabel("Close")
plt.legend()
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

[cell 72 markdown]
### Quick demo: next-week forecast function for any ticker

This helper takes the last `SEQ_LEN` weeks of OHLCV for a ticker and returns a single next-week close forecast.

[cell 73 code]
def forecast_next_week_close(ticker: str):
    if ticker not in stock_to_id:
        raise ValueError(f"Unknown ticker: {ticker}")

    g = df_v1[df_v1["stock"] == ticker].sort_values("date").reset_index(drop=True)
    if len(g) <= SEQ_LEN:
        raise ValueError("Not enough history for the selected SEQ_LEN")

    x_last = g[FEATURE_COLS].astype(float).values[-SEQ_LEN:]
    x_last_scaled = scaler_X.transform(x_last)
    x_t = torch.tensor(x_last_scaled.reshape(1, SEQ_LEN, len(FEATURE_COLS)), dtype=torch.float32).to(device)
    sid_t = torch.tensor([stock_to_id[ticker]], dtype=torch.long).to(device)

    with torch.no_grad():
        pred_scaled = model(x_t, sid_t).cpu().numpy()
    pred = float(scaler_y.inverse_transform(pred_scaled)[0, 0])
    return pred

example_ticker = stocks[0]
print(example_ticker, "-> forecast next-week close:", forecast_next_week_close(example_ticker))

[cell 74 markdown]
## LSTM demos

- Dow Jones multi-horizon forecasting with an LSTM (direct vs recursive decoding)
- Walk-forward evaluation on one ticker
- EEG eye-state modeling with an LSTM (many-to-one and many-to-many)
- Basic interpretability via input-gradient saliency

[cell 75 markdown]
### 0) Setup

[cell 76 code]
import os
import math
import random

import numpy as np
import pandas as pd

import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader, TensorDataset

from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error, accuracy_score, f1_score, roc_auc_score, confusion_matrix

import matplotlib.pyplot as plt

SEED = 42
random.seed(SEED)
np.random.seed(SEED)
torch.manual_seed(SEED)

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
device

[cell 77 markdown]
### 1) Dow Jones: multi-horizon forecasting with a vanilla LSTM (direct prediction)

[cell 78 code]
# Data load + cleanup (self-contained; does not depend on earlier cells)
DJ_PATH = "dow_jones_index.data"
if not os.path.exists(DJ_PATH):
    DJ_PATH = "/mnt/data/dow_jones_index.data"

dj = pd.read_csv(DJ_PATH)

def _money_to_float(x):
    if isinstance(x, str):
        return float(x.replace("$", "").replace(",", ""))
    return float(x)

for c in ["open", "high", "low", "close", "next_weeks_open", "next_weeks_close"]:
    dj[c] = dj[c].apply(_money_to_float)

dj["date"] = pd.to_datetime(dj["date"])
dj = dj.sort_values(["stock", "date"]).reset_index(drop=True)

# Features that do not leak future information (exclude next_weeks_* and percent_change_next_weeks_price)
feature_cols = [
    "open", "high", "low", "close", "volume",
    "percent_change_price",
    "percent_change_volume_over_last_wk",
    "previous_weeks_volume",
    "days_to_next_dividend",
    "percent_return_next_dividend",
]

dj = dj.dropna(subset=feature_cols).reset_index(drop=True)

dj[["stock", "date"]].head(), dj.shape

[cell 79 code]
# Train/val/test split per ticker (chronological), and global scaling fit on train only
train_frac = 0.80
val_frac_within_train = 0.10  # last part of train window

train_parts = []
for stock, g in dj.groupby("stock", sort=False):
    g = g.sort_values("date")
    split = int(len(g) * train_frac)
    train_parts.append(g.iloc[:split])

dj_train_rows = pd.concat(train_parts, axis=0)
scaler = StandardScaler().fit(dj_train_rows[feature_cols])

dj_scaled = dj.copy()
dj_scaled[feature_cols] = scaler.transform(dj_scaled[feature_cols])

# Sequence builder: predict future close returns relative to the last close in the input window
LOOKBACK = 8   # weeks
HORIZON = 4    # weeks ahead (predict 1..HORIZON)

stock_list = sorted(dj_scaled["stock"].unique())
stock_to_id = {s: i for i, s in enumerate(stock_list)}

X_tr, y_tr, sid_tr = [], [], []
X_va, y_va, sid_va = [], [], []
X_te, y_te, sid_te = [], [], []

# Keep metadata for plotting
meta_te = []  # (stock, end_date, base_close, future_dates, future_close)

for stock, g in dj_scaled.groupby("stock", sort=False):
    g = g.sort_values("date").reset_index(drop=True)
    split = int(len(g) * train_frac)
    val_cut = int(split * (1.0 - val_frac_within_train))

    X = g[feature_cols].values.astype(np.float32)
    close = g["close"].values.astype(np.float32)
    dates = g["date"].values
    sid = stock_to_id[stock]

    for t in range(LOOKBACK - 1, len(g) - HORIZON):
        x_seq = X[t - LOOKBACK + 1 : t + 1]
        base_close = close[t]
        future_close = close[t + 1 : t + HORIZON + 1]
        y_seq = (future_close / base_close) - 1.0  # future returns relative to base_close

        target_last_idx = t + HORIZON
        if target_last_idx < val_cut:
            X_tr.append(x_seq); y_tr.append(y_seq); sid_tr.append(sid)
        elif target_last_idx < split:
            X_va.append(x_seq); y_va.append(y_seq); sid_va.append(sid)
        else:
            X_te.append(x_seq); y_te.append(y_seq); sid_te.append(sid)
            meta_te.append((stock, pd.to_datetime(dates[t]), float(base_close),
                            pd.to_datetime(dates[t+1:t+HORIZON+1]), future_close.astype(float)))

X_tr = np.stack(X_tr); y_tr = np.stack(y_tr); sid_tr = np.array(sid_tr)
X_va = np.stack(X_va); y_va = np.stack(y_va); sid_va = np.array(sid_va)
X_te = np.stack(X_te); y_te = np.stack(y_te); sid_te = np.array(sid_te)

X_tr.shape, y_tr.shape, X_va.shape, X_te.shape

[cell 80 code]
class SeqDataset(Dataset):
    def __init__(self, X: np.ndarray, y: np.ndarray):
        self.X = torch.tensor(X, dtype=torch.float32)
        self.y = torch.tensor(y, dtype=torch.float32)

    def __len__(self):
        return self.X.shape[0]

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

train_loader = DataLoader(SeqDataset(X_tr, y_tr), batch_size=128, shuffle=True)
val_loader   = DataLoader(SeqDataset(X_va, y_va), batch_size=256, shuffle=False)
test_loader  = DataLoader(SeqDataset(X_te, y_te), batch_size=256, shuffle=False)

len(train_loader), len(val_loader), len(test_loader)

[cell 81 code]
class LSTMMultiHorizon(nn.Module):
    def __init__(self, n_features: int, hidden_size: int, num_layers: int, horizon: int, dropout: float = 0.2):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=(dropout if num_layers > 1 else 0.0),
        )
        self.head = nn.Linear(hidden_size, horizon)

    def forward(self, x):
        _, (h, _) = self.lstm(x)
        return self.head(h[-1])

def train_regression(model, train_loader, val_loader, epochs=30, lr=1e-3, clip_norm=1.0, patience=5):
    model = model.to(device)
    opt = torch.optim.Adam(model.parameters(), lr=lr)
    sched = torch.optim.lr_scheduler.ReduceLROnPlateau(opt, mode="min", factor=0.5, patience=2)
    loss_fn = nn.MSELoss()

    best_val = float("inf")
    best_state = None
    wait = 0
    history = {"train_loss": [], "val_loss": []}

    for epoch in range(1, epochs + 1):
        model.train()
        tr_losses = []
        for xb, yb in train_loader:
            xb = xb.to(device)
            yb = yb.to(device)

            pred = model(xb)
            loss = loss_fn(pred, yb)

            opt.zero_grad()
            loss.backward()
            if clip_norm is not None:
                nn.utils.clip_grad_norm_(model.parameters(), clip_norm)
            opt.step()

            tr_losses.append(loss.item())

        model.eval()
        va_losses = []
        with torch.no_grad():
            for xb, yb in val_loader:
                xb = xb.to(device)
                yb = yb.to(device)
                va_losses.append(loss_fn(model(xb), yb).item())

        tr = float(np.mean(tr_losses)) if tr_losses else float("nan")
        va = float(np.mean(va_losses)) if va_losses else float("nan")
        history["train_loss"].append(tr)
        history["val_loss"].append(va)

        sched.step(va)
        print(f"Epoch {epoch:02d} | train MSE {tr:.6f} | val MSE {va:.6f} | lr {opt.param_groups[0]['lr']:.2e}")

        if va < best_val - 1e-6:
            best_val = va
            best_state = {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}
            wait = 0
        else:
            wait += 1
            if wait >= patience:
                print("Early stopping.")
                break

    if best_state is not None:
        model.load_state_dict(best_state)
    return model, history

model_mh = LSTMMultiHorizon(n_features=len(feature_cols), hidden_size=96, num_layers=2, horizon=HORIZON, dropout=0.2)
model_mh, hist_mh = train_regression(model_mh, train_loader, val_loader, epochs=25, lr=1e-3, clip_norm=1.0, patience=5)

[cell 82 code]
plt.figure()
plt.plot(hist_mh["train_loss"], label="train")
plt.plot(hist_mh["val_loss"], label="val")
plt.xlabel("epoch")
plt.ylabel("MSE")
plt.title("Dow Jones LSTM multi-horizon training curve")
plt.legend()
plt.show()

[cell 83 code]
def eval_multi_horizon(model, loader):
    model.eval()
    ys, ps = [], []
    with torch.no_grad():
        for xb, yb in loader:
            xb = xb.to(device)
            pred = model(xb).cpu().numpy()
            ys.append(yb.numpy())
            ps.append(pred)
    y = np.concatenate(ys, axis=0)
    p = np.concatenate(ps, axis=0)
    return y, p

y_true, y_pred = eval_multi_horizon(model_mh.to(device), test_loader)

h_metrics = []
for h in range(HORIZON):
    mae = mean_absolute_error(y_true[:, h], y_pred[:, h])
    rmse = math.sqrt(mean_squared_error(y_true[:, h], y_pred[:, h]))
    h_metrics.append((h+1, mae, rmse))

pd.DataFrame(h_metrics, columns=["horizon_weeks", "MAE_return", "RMSE_return"])

[cell 84 code]
# Plot one ticker's test predictions: reconstruct price paths from predicted returns
plot_ticker = stock_list[0]

idxs = [i for i, (s, *_rest) in enumerate(meta_te) if s == plot_ticker]
if len(idxs) == 0:
    plot_ticker = meta_te[0][0]
    idxs = [i for i, (s, *_rest) in enumerate(meta_te) if s == plot_ticker]

# Take a handful of samples for plotting
idxs = idxs[:15]

rows = []
for j in idxs:
    stock, end_date, base_close, future_dates, future_close = meta_te[j]
    pred_returns = y_pred[j]
    pred_close = base_close * (1.0 + pred_returns)

    rows.append({
        "end_date": end_date,
        "base_close": base_close,
        "true_close_t+1": float(future_close[0]),
        "pred_close_t+1": float(pred_close[0]),
    })

pd.DataFrame(rows).head()

[cell 85 code]
# Visualize multi-step paths for a single sample
sample_idx = idxs[0]
stock, end_date, base_close, future_dates, future_close = meta_te[sample_idx]
pred_returns = y_pred[sample_idx]
pred_close = base_close * (1.0 + pred_returns)

x_dates = [end_date] + list(future_dates)
true_path = [base_close] + list(future_close)
pred_path  = [base_close] + list(pred_close)

plt.figure()
plt.plot(x_dates, true_path, marker="o", label="true")
plt.plot(x_dates, pred_path, marker="o", label="pred")
plt.title(f"LSTM direct multi-horizon forecast path ({stock})")
plt.xlabel("date")
plt.ylabel("close")
plt.legend()
plt.xticks(rotation=30)
plt.show()

[cell 86 markdown]
### 2) Direct vs recursive multi-step forecasting with the same vanilla LSTM (different decoding)

[cell 87 code]
# One-step model trained to predict only the 1-week-ahead return.
# Then we decode recursively to reach a multi-week horizon.

class LSTMOneStep(nn.Module):
    def __init__(self, n_features: int, hidden_size: int, num_layers: int, dropout: float = 0.2):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=(dropout if num_layers > 1 else 0.0),
        )
        self.head = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (h, _) = self.lstm(x)
        return self.head(h[-1]).squeeze(1)

# Prepare one-step targets (h=1) from the already-built arrays
y_tr_1 = y_tr[:, 0:1]
y_va_1 = y_va[:, 0:1]
y_te_1 = y_te[:, 0:1]

train1 = DataLoader(SeqDataset(X_tr, y_tr_1), batch_size=128, shuffle=True)
val1   = DataLoader(SeqDataset(X_va, y_va_1), batch_size=256, shuffle=False)
test1  = DataLoader(SeqDataset(X_te, y_te_1), batch_size=256, shuffle=False)

model_1 = LSTMOneStep(n_features=len(feature_cols), hidden_size=96, num_layers=2, dropout=0.2)
model_1, hist_1 = train_regression(model_1, train1, val1, epochs=15, lr=1e-3, clip_norm=1.0, patience=4)

[cell 88 code]
# Recursive decoding helper:
# We keep the input window fixed except for the last 'close' feature.
# This is a demo-friendly approximation; the goal is to show the decoding idea, not perfect market simulation.

close_ix = feature_cols.index("close")

def recursive_forecast_returns(model, x_seq_scaled: np.ndarray, base_close: float, steps: int = 4):
    model.eval()
    x = torch.tensor(x_seq_scaled[None, :, :], dtype=torch.float32).to(device)
    cur_close = float(base_close)
    preds = []

    for _ in range(steps):
        with torch.no_grad():
            r = float(model(x).cpu().numpy().ravel()[0])
        preds.append(r)
        cur_close = cur_close * (1.0 + r)

        # Update only the last timestep's 'close' feature to reflect the predicted close (scaled).
        x_np = x.cpu().numpy()
        tmp_scaled = (cur_close - scaler.mean_[close_ix]) / scaler.scale_[close_ix]
        x_np[0, -1, close_ix] = tmp_scaled
        x = torch.tensor(x_np, dtype=torch.float32).to(device)

    return np.array(preds, dtype=np.float32)

# Compare direct multi-horizon vs recursive for a few test samples
compare_n = min(200, len(X_te))
direct = y_pred[:compare_n]                      # from multi-horizon model
recur = np.stack([recursive_forecast_returns(model_1.to(device), X_te[i], meta_te[i][2], HORIZON)
                  for i in range(compare_n)], axis=0)
truth = y_true[:compare_n]

rows = []
for h in range(HORIZON):
    rows.append({
        "horizon": h+1,
        "MAE_direct": mean_absolute_error(truth[:, h], direct[:, h]),
        "MAE_recursive": mean_absolute_error(truth[:, h], recur[:, h]),
    })

pd.DataFrame(rows)

[cell 89 markdown]
### 3) Walk-forward evaluation on one ticker (expanding window retrain)

[cell 90 code]
def walk_forward_one_ticker(dj_scaled: pd.DataFrame, ticker: str, lookback: int, horizon: int, epochs_per_step: int = 8):
    g = dj_scaled[dj_scaled["stock"] == ticker].sort_values("date").reset_index(drop=True)
    X = g[feature_cols].values.astype(np.float32)
    close = g["close"].values.astype(np.float32)
    dates = g["date"].values

    start_train_end = max(lookback + 5, int(len(g) * 0.6))
    preds, trues, pred_dates = [], [], []

    for train_end in range(start_train_end, len(g) - horizon):
        # Build train sequences up to train_end (targets <= train_end)
        X_tr_w, y_tr_w = [], []
        for t in range(lookback - 1, train_end - horizon + 1):
            x_seq = X[t - lookback + 1 : t + 1]
            base_close = close[t]
            future_close = close[t + 1 : t + horizon + 1]
            y_seq = (future_close / base_close) - 1.0
            X_tr_w.append(x_seq); y_tr_w.append(y_seq)

        if len(X_tr_w) < 20:
            continue

        X_tr_w = np.stack(X_tr_w); y_tr_w = np.stack(y_tr_w)
        loader = DataLoader(SeqDataset(X_tr_w, y_tr_w), batch_size=64, shuffle=True)

        model = LSTMMultiHorizon(n_features=len(feature_cols), hidden_size=64, num_layers=1, horizon=horizon, dropout=0.0).to(device)
        opt = torch.optim.Adam(model.parameters(), lr=1e-3)
        loss_fn = nn.MSELoss()

        for _ in range(epochs_per_step):
            model.train()
            for xb, yb in loader:
                xb = xb.to(device); yb = yb.to(device)
                pred = model(xb)
                loss = loss_fn(pred, yb)
                opt.zero_grad()
                loss.backward()
                nn.utils.clip_grad_norm_(model.parameters(), 1.0)
                opt.step()

        # Forecast from the last available window at train_end
        t = train_end - 1
        x_seq = X[t - lookback + 1 : t + 1]
        base_close = float(close[t])
        true_future = close[t + 1 : t + horizon + 1]
        true_returns = (true_future / base_close) - 1.0

        with torch.no_grad():
            pred_returns = model(torch.tensor(x_seq[None, :, :], dtype=torch.float32).to(device)).cpu().numpy().ravel()

        preds.append(pred_returns)
        trues.append(true_returns)
        pred_dates.append(pd.to_datetime(dates[t]))

    return np.array(preds), np.array(trues), pd.to_datetime(pred_dates)

wf_ticker = stock_list[0]
wf_pred, wf_true, wf_dates = walk_forward_one_ticker(dj_scaled, wf_ticker, LOOKBACK, HORIZON, epochs_per_step=6)

wf_pred.shape, wf_true.shape, wf_ticker

[cell 91 code]
# Plot walk-forward 1-step-ahead returns (horizon=1 slice)
if len(wf_dates) > 0:
    plt.figure()
    plt.plot(wf_dates, wf_true[:, 0], label="true")
    plt.plot(wf_dates, wf_pred[:, 0], label="pred")
    plt.title(f"Walk-forward 1-week return prediction ({wf_ticker})")
    plt.xlabel("date")
    plt.ylabel("return")
    plt.legend()
    plt.xticks(rotation=30)
    plt.show()

    print("Walk-forward MAE (h=1):", mean_absolute_error(wf_true[:, 0], wf_pred[:, 0]))
else:
    print("Walk-forward run produced no points (ticker too short after cleaning).")

[cell 92 markdown]
### 4) EEG: LSTM sequence classifier (many-to-one) using stride-based windows

[cell 93 code]
from scipy.io import arff

EEG_PATH = "EEG Eye State.arff"
if not os.path.exists(EEG_PATH):
    EEG_PATH = "/mnt/data/EEG Eye State.arff"

eeg_raw, eeg_meta = arff.loadarff(EEG_PATH)
eeg = pd.DataFrame(eeg_raw)

# eyeDetection comes as bytes b'0'/'1'
eeg["eyeDetection"] = eeg["eyeDetection"].apply(lambda x: int(x.decode("utf-8")))

eeg_feature_cols = [c for c in eeg.columns if c != "eyeDetection"]
X = eeg[eeg_feature_cols].values.astype(np.float32)
y = eeg["eyeDetection"].values.astype(np.int64)

# Chronological split (EEG rows are time-ordered)
split = int(len(eeg) * 0.8)
X_train, X_test = X[:split], X[split:]
y_train, y_test = y[:split], y[split:]

sc_eeg = StandardScaler().fit(X_train)
X_train = sc_eeg.transform(X_train).astype(np.float32)
X_test  = sc_eeg.transform(X_test).astype(np.float32)

SEQ_LEN = 64
STRIDE = 8

def make_windows(X_arr, y_arr, seq_len=64, stride=8):
    Xs, ys = [], []
    for start in range(0, len(X_arr) - seq_len + 1, stride):
        end = start + seq_len
        Xs.append(X_arr[start:end])
        ys.append(y_arr[end - 1])  # label at the end of the window
    return np.stack(Xs), np.array(ys)

Xtr_w, ytr_w = make_windows(X_train, y_train, SEQ_LEN, STRIDE)
Xte_w, yte_w = make_windows(X_test, y_test, SEQ_LEN, STRIDE)

Xtr_w.shape, ytr_w.mean(), Xte_w.shape, yte_w.mean()

[cell 94 code]
class LSTMClassifier(nn.Module):
    def __init__(self, n_features: int, hidden_size: int = 64, num_layers: int = 1, dropout: float = 0.0):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
            dropout=(dropout if num_layers > 1 else 0.0),
        )
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        _, (h, _) = self.lstm(x)
        return self.fc(h[-1]).squeeze(1)

def train_classifier(model, Xtr, ytr, Xva, yva, epochs=15, lr=1e-3, batch_size=128, clip_norm=1.0, patience=4):
    model = model.to(device)
    loss_fn = nn.BCEWithLogitsLoss()
    opt = torch.optim.Adam(model.parameters(), lr=lr)

    train_loader = DataLoader(SeqDataset(Xtr, ytr.reshape(-1, 1).astype(np.float32)), batch_size=batch_size, shuffle=True)
    val_loader   = DataLoader(SeqDataset(Xva, yva.reshape(-1, 1).astype(np.float32)), batch_size=256, shuffle=False)

    best_val = float("inf")
    best_state = None
    wait = 0

    for epoch in range(1, epochs + 1):
        model.train()
        for xb, yb in train_loader:
            xb = xb.to(device)
            yb = yb.to(device).squeeze(1)
            logits = model(xb)
            loss = loss_fn(logits, yb)

            opt.zero_grad()
            loss.backward()
            if clip_norm is not None:
                nn.utils.clip_grad_norm_(model.parameters(), clip_norm)
            opt.step()

        model.eval()
        losses = []
        with torch.no_grad():
            for xb, yb in val_loader:
                xb = xb.to(device)
                yb = yb.to(device).squeeze(1)
                losses.append(loss_fn(model(xb), yb).item())
        val_loss = float(np.mean(losses))

        print(f"Epoch {epoch:02d} | val BCE {val_loss:.5f}")
        if val_loss < best_val - 1e-6:
            best_val = val_loss
            best_state = {k: v.detach().cpu().clone() for k, v in model.state_dict().items()}
            wait = 0
        else:
            wait += 1
            if wait >= patience:
                print("Early stopping.")
                break

    if best_state is not None:
        model.load_state_dict(best_state)
    return model

# Validation split from the end of train windows
val_cut = int(len(Xtr_w) * 0.9)
Xtr2, ytr2 = Xtr_w[:val_cut], ytr_w[:val_cut]
Xva2, yva2 = Xtr_w[val_cut:], ytr_w[val_cut:]

clf = LSTMClassifier(n_features=len(eeg_feature_cols), hidden_size=64, num_layers=1, dropout=0.0)
clf = train_classifier(clf, Xtr2, ytr2, Xva2, yva2, epochs=20, lr=1e-3, batch_size=128, clip_norm=1.0, patience=4)

[cell 95 code]
# EEG evaluation
clf.eval()
with torch.no_grad():
    logits = clf(torch.tensor(Xte_w, dtype=torch.float32).to(device)).cpu().numpy()
proba = 1.0 / (1.0 + np.exp(-logits))
pred = (proba >= 0.5).astype(int)

acc = accuracy_score(yte_w, pred)
f1 = f1_score(yte_w, pred)
auc = roc_auc_score(yte_w, proba)

cm = confusion_matrix(yte_w, pred)

print({"accuracy": acc, "f1": f1, "roc_auc": auc})
cm

[cell 96 code]
# Confusion matrix plot
plt.figure()
plt.imshow(cm)
plt.title("EEG LSTM confusion matrix")
plt.xlabel("pred")
plt.ylabel("true")
plt.xticks([0,1])
plt.yticks([0,1])
for (i, j), v in np.ndenumerate(cm):
    plt.text(j, i, str(v), ha="center", va="center")
plt.show()

[cell 97 markdown]
### 5) EEG: many-to-many sequence labeling with a vanilla LSTM

[cell 98 code]
# Many-to-many: predict eye state for each timestep in the window.
# This shows a different LSTM usage mode (sequence labeling) while staying with vanilla LSTM.

def make_windows_seq_labels(X_arr, y_arr, seq_len=64, stride=32):
    Xs, Ys = [], []
    for start in range(0, len(X_arr) - seq_len + 1, stride):
        end = start + seq_len
        Xs.append(X_arr[start:end])
        Ys.append(y_arr[start:end])
    return np.stack(Xs), np.stack(Ys)

Xtr_seq, ytr_seq = make_windows_seq_labels(X_train, y_train, seq_len=SEQ_LEN, stride=32)
Xte_seq, yte_seq = make_windows_seq_labels(X_test, y_test, seq_len=SEQ_LEN, stride=32)

Xtr_seq.shape, ytr_seq.shape

[cell 99 code]
class LSTMTagger(nn.Module):
    def __init__(self, n_features: int, hidden_size: int = 64, num_layers: int = 1):
        super().__init__()
        self.lstm = nn.LSTM(
            input_size=n_features,
            hidden_size=hidden_size,
            num_layers=num_layers,
            batch_first=True,
        )
        self.fc = nn.Linear(hidden_size, 1)

    def forward(self, x):
        out, _ = self.lstm(x)              # (B, T, H)
        logits = self.fc(out).squeeze(-1)  # (B, T)
        return logits

def train_tagger(model, Xtr, ytr, epochs=8, lr=1e-3, batch_size=128):
    model = model.to(device)
    opt = torch.optim.Adam(model.parameters(), lr=lr)
    loss_fn = nn.BCEWithLogitsLoss()

    ds = TensorDataset(
        torch.tensor(Xtr, dtype=torch.float32),
        torch.tensor(ytr, dtype=torch.float32),
    )
    loader = DataLoader(ds, batch_size=batch_size, shuffle=True)

    for epoch in range(1, epochs + 1):
        model.train()
        losses = []
        for xb, yb in loader:
            xb = xb.to(device)
            yb = yb.to(device)

            logits = model(xb)
            loss = loss_fn(logits, yb)

            opt.zero_grad()
            loss.backward()
            nn.utils.clip_grad_norm_(model.parameters(), 1.0)
            opt.step()

            losses.append(loss.item())

        print(f"Epoch {epoch:02d} | train BCE {float(np.mean(losses)):.5f}")
    return model

tagger = LSTMTagger(n_features=len(eeg_feature_cols), hidden_size=64, num_layers=1)
tagger = train_tagger(tagger, Xtr_seq, ytr_seq, epochs=8, lr=1e-3, batch_size=128)

[cell 100 code]
# Evaluate tagger on test sequences
tagger.eval()
with torch.no_grad():
    logits = tagger(torch.tensor(Xte_seq, dtype=torch.float32).to(device)).cpu().numpy()
proba = 1.0 / (1.0 + np.exp(-logits))
pred = (proba >= 0.5).astype(int)

# Flatten across time for a quick summary metric
flat_true = yte_seq.reshape(-1)
flat_pred = pred.reshape(-1)
flat_proba = proba.reshape(-1)

print({
    "accuracy": accuracy_score(flat_true, flat_pred),
    "f1": f1_score(flat_true, flat_pred),
    "roc_auc": roc_auc_score(flat_true, flat_proba),
})

[cell 101 code]
# Visualize one labeled window (true vs predicted probabilities across time)
win = 0
t = np.arange(SEQ_LEN)

plt.figure()
plt.plot(t, yte_seq[win], label="true")
plt.plot(t, proba[win], label="pred_proba")
plt.title("EEG LSTM sequence labeling (one window)")
plt.xlabel("timestep")
plt.ylabel("eye_state / probability")
plt.legend()
plt.show()

[cell 102 markdown]
### 6) LSTM interpretability: input-gradient saliency (Dow Jones)

[cell 103 code]
# Saliency: gradient of a chosen output w.r.t. input features over timesteps.
# We compute this for the multi-horizon forecaster and visualize as a heatmap.

model_mh.eval()

# Pick one test sample
i = 0
x_np = X_te[i].copy()
x = torch.tensor(x_np[None, :, :], dtype=torch.float32, requires_grad=True).to(device)

# Choose horizon step 1 output
out = model_mh(x)[0, 0]
out.backward()

sal = x.grad.detach().cpu().numpy()[0]  # (T, F)
sal = np.abs(sal)

plt.figure()
plt.imshow(sal, aspect="auto")
plt.title("Saliency |d(output_h=1)/d(inputs)|")
plt.xlabel("feature")
plt.ylabel("timestep")
plt.xticks(np.arange(len(feature_cols)), feature_cols, rotation=90)
plt.colorbar()
plt.tight_layout()
plt.show()

[cell 104 markdown]
### Summary

- For forecasting, show the training curve, horizon-wise errors, and the multi-step path plot.
- For decoding, compare direct multi-horizon vs recursive decoding (table of MAE by horizon).
- For evaluation, use walk-forward to explain why chronological testing matters.
- For EEG, contrast many-to-one vs many-to-many LSTM usage.
- For interpretability, use saliency heatmaps to explain what timesteps/features the LSTM is sensitive to.