## Introduction

In this tutorial, we will explore several machine learning classification models to predict the direction of stock movement for the **Pakistan Stock Exchange (PSX)** using historical data. We will use models such as **Logistic Regression**, **Linear Discriminant Analysis (LDA)**, **Quadratic Discriminant Analysis (QDA)**, **K-Nearest Neighbors (KNN)**, and **Decision Trees** (CART).

The goal is to understand some machine learning techniques and how they can be applied to real-world financial data. We will use historical stock data from the PSX. We can understand that predicting movements in stock prices is a challenging task due to the complex nature of financial markets. This is a simplified example to demonstrate how machine learning models can be used in real financial applications. We have created confusion matrix and calculated accuracy, precision, recall, and F1 score to evaluate the performance of each model. We
predict whether the stock price will go "Up" or "Down" based on historical changes in stock price and trading volume. 


## 1. Loading and Preparing the Data

The first step is to load the PSX data, clean it, and perform necessary transformations for feature engineering.

```{r}
#| label: load-libraries
library(caret)      # For model evaluation and cross-validation
library(rpart)      # For Decision Trees
library(tidyverse)  # For data manipulation and visualization
library(Metrics)    # For additional evaluation metrics
library(janitor)    # For cleaning column names
library(MASS)
library(class)    # For KNN
```


## Load data

We have psx data in csv format. We will load the data and inspect its structure.


```{r}
#| label: load-data
psx <- read_csv("data/psx.csv")

# Inspect the data structure (rows, columns, and types)
dim(psx)
glimpse(psx)
```

## Data cleaning and feature engineering

We will clean the data, convert columns to appropriate types, create new features, and handle missing values.

```{r}
#| label: data-cleaning

# Convert the `Date` column to a proper Date format
psx <- psx %>% 
  mutate(date = as.Date(Date, format = "%m/%d/%Y")) |> dplyr::select(-Date)
# Clean column names to make them consistent and readable
psx <- psx %>% clean_names()

# Convert `change_percent` to numeric and create a new column `direction` for the stock movement (Up/Down)
psx <- psx %>% 
  mutate(
    change_percent = as.numeric(str_remove(change_percent, "%")),
    direction = if_else(change_percent > 0, "Up", "Down")
  )

# Arrange data by ascending date to maintain temporal order
psx <- psx %>% arrange(date)

# Create lag variables for `change_percent` (lags of 1 to 5 days)
psx <- psx %>% 
  mutate(
    lag1 = lag(change_percent, 1),
    lag2 = lag(change_percent, 2),
    lag3 = lag(change_percent, 3),
    lag4 = lag(change_percent, 4),
    lag5 = lag(change_percent, 5)
  )

# Convert `vol` (volume) to numeric, handling the M and K suffixes
psx <- psx %>% mutate(vol = as.numeric(str_replace_all(vol, c("M" = "e6", "K" = "e3"))))

# Drop rows with missing values (NA)
psx_clean <- psx %>% drop_na()

# Convert `direction` to numeric for binary classification (1 for "Up", 0 for "Down")
psx_clean <- psx_clean %>%
  mutate(direction = if_else(direction == "Up", 1, 0))
```

There are `r nrow(psx_clean)` rows in the cleaned dataset after data cleaning and feature engineering.


**Explanation:** 

1. We loaded the dataset and cleaned column names. 

2. We transformed the `change_percent` column from character to numeric. 

3. We created a `direction` column to classify stock movements as "Up" or "Down". 

4. Lag features are created to provide past data (1-day lag to 5-day lag) as predictors for stock movement. 

5. We converted the volume (`vol`) column to numeric format for use in the model. 

6. Rows with missing values were dropped, and the target variable `direction` was converted to binary numeric values for classification.



## 2. Logistic Regression (GLM)

Logistic regression will model the probability that the stock will go "Up" based on the lag features and volume.

```{r}
#| evaluate: false
# Fit logistic regression model
# glm_fit <- glm(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, 
#                data = psx_clean, family = binomial)
# 
# # Model summary
# summary(glm_fit)
# 
# # Add predicted probabilities to the dataset
# psx_clean <- psx_clean %>%
#   mutate(
#     glm_probs = predict(glm_fit, type = "response"),
#     glm_pred = if_else(glm_probs > 0.5, 1, 0)
#   )
# 
# # Evaluate logistic regression performance using confusion matrix and accuracy
# confusion_matrix <- table(glm_pred = psx_clean$glm_pred, actual = psx_clean$direction)
# accuracy <- mean(psx_clean$glm_pred == psx_clean$direction)

# Display results
#print(confusion_matrix)
#print(paste("Overall Accuracy:", round(accuracy, 3)))
```

**Explanation:**
1. The logistic regression model is fit with the lag features and volume as predictors.
2. Predictions are made for each observation using `predict()`, and the predicted probabilities are added to the dataset.
3. We then classify the probability as "Up" or "Down" based on a threshold of 0.5.
4. The performance is evaluated by comparing predicted values against actual values using a confusion matrix and accuracy.

---

## 3. Train-Test Split

We split the data into **training** and **testing** datasets to evaluate the performance of our models.

```{r}
# Train-test split based on the date
train_data <- psx_clean %>%
  filter(date < as.Date("2024-01-01"))

test_data <- psx_clean %>%
  filter(date >= as.Date("2024-01-01"))
# 
# # Ensure direction is numeric (0 for Down, 1 for Up)
# train_data <- train_data %>%
#   mutate(direction = as.factor(direction))
# 
# test_data <- test_data %>%
#   mutate(direction = as.factor(direction))
# 
# # Ensure consistent factor levels for direction in both training and test data
# train_data$direction <- factor(train_data$direction, levels = c(0, 1))  # 0 = Down, 1 = Up
# test_data$direction <- factor(test_data$direction, levels = c(0, 1))

```



## Linear Discriminant Analysis (LDA)


```{r}
#| label: lda-model

lda_fit <- lda(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol,
               data = train_data)

# Predict on test data

lda_pred <- predict(lda_fit, newdata = test_data)

# Evaluate LDA model

lda_accuracy <- mean(lda_pred$class == test_data$direction)
```

**Explanation:**
1. We split the data based on the date (all data before January 1, 2024 for training, and the rest for testing).
2. The target variable `direction` is converted into a factor to comply with classification model requirements.

---

## 4. Linear Discriminant Analysis (LDA)

LDA is another classification technique that assumes normality in the predictor variables and tries to separate the classes (Up/Down) by maximizing the ratio of between-class variance to within-class variance.



Error in Ops.factor(lda_pred, test_data$direction) : 
  level sets of factors are different
> 



**Explanation:**
1. We fit the LDA model using the training data.
2. Predictions are made on the test set using the fitted model.
3. The model accuracy is evaluated by comparing predicted and actual values.

---

## 5. Quadratic Discriminant Analysis (QDA)

QDA is similar to LDA but allows for each class to have its own covariance matrix, making it more flexible when the data distribution is not the same across classes.

```{r}
#| label: qda-model
# Fit QDA model
qda_fit <- qda(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, data = train_data)

# Predict on test data
qda_pred <- predict(qda_fit, newdata = test_data)$class

# Evaluate QDA model performance
qda_accuracy <- mean(qda_pred == test_data$direction)

# Display QDA results
print(paste("QDA Accuracy:", round(qda_accuracy, 3)))
```

**Explanation:**
1. We fit the QDA model using the training data.
2. Predictions are made on the test set, and the model’s accuracy is evaluated.

---

## 6. K-Nearest Neighbors (KNN)

KNN is a non-parametric method that assigns the class of an observation based on the majority class of its k nearest neighbors in the feature space.

```{r}
# Fit KNN model

knn_fit <- knn(train = train_data[, c("lag1", "lag2", "lag3", "lag4", "lag5", "vol")],
               test = test_data[, c("lag1", "lag2", "lag3", "lag4", "lag5", "vol")],
               cl = train_data$direction,
               k = 5)

# Evaluate KNN model

knn_accuracy <- mean(knn_fit == test_data$direction)
```

**Explanation:**
1. Data is normalized to ensure all features have the same scale for KNN.
2. We fit the KNN model with `k = 5` and evaluate its accuracy.

---

## 7. Decision Trees (CART)

The Decision Tree algorithm creates a tree-like model of decisions and their possible consequences.

```{r}
# Fit Decision Tree model (CART)
dt_fit <- rpart(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol,
                 data = train_data,
                 method = "class")

# Predict on test data

dt_pred <- predict(dt_fit, newdata = test_data, type = "class")


# Evaluate Decision Tree model

dt_accuracy <- mean(dt_pred == test_data$direction)

```

**Explanation:**
1. We fit the Decision Tree (CART) model using the training data.
2. Predictions are made on the test set, and the accuracy is calculated.

## Logistic Regression

```{r}
# First ensure direction is properly encoded as numeric binary in psx_clean
psx_clean <- psx_clean %>%
  mutate(direction = as.numeric(direction))  # Should be 0 or 1

# Logistic regression model
glm_fit <- glm(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol,
               data = psx_clean, family = binomial)

# Model summary
summary(glm_fit)

# Predict probabilities and classify directions
psx_clean <- psx_clean %>%
  mutate(
    glm_probs = predict(glm_fit, type = "response"),
    glm_pred = as.numeric(glm_probs > 0.5)  # Convert to numeric 0/1
  )

# Evaluate model accuracy
confusion_matrix <- table(Predicted = psx_clean$glm_pred, Actual = psx_clean$direction)
accuracy <- mean(psx_clean$glm_pred == psx_clean$direction)

# Split data into train and test sets (based on date < 2024)
train <- psx_clean %>%
  filter(date < as.Date("2024-01-01"))

test <- psx_clean %>%
  filter(date >= as.Date("2024-01-01"))

# Train logistic regression on train data
glm_train <- glm(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol,
                 data = train, family = binomial)

# Predict on test data
test <- test %>%
  mutate(
    glm_probs = predict(glm_train, newdata = ., type = "response"),
    glm_pred = as.numeric(glm_probs > 0.5)  # Convert to numeric 0/1
  )

# Evaluate test accuracy
confusion_matrix_test <- table(Predicted = test$glm_pred, Actual = test$direction)
test_accuracy <- mean(test$glm_pred == test$direction)

# Create detailed evaluation metrics
train_metrics <- list(
  accuracy = accuracy,
  sensitivity = confusion_matrix[2,2] / sum(confusion_matrix[,2]),
  specificity = confusion_matrix[1,1] / sum(confusion_matrix[,1]),
  confusion_matrix = confusion_matrix
)

test_metrics <- list(
  accuracy = test_accuracy,
  sensitivity = confusion_matrix_test[2,2] / sum(confusion_matrix_test[,2]),
  specificity = confusion_matrix_test[1,1] / sum(confusion_matrix_test[,1]),
  confusion_matrix = confusion_matrix_test
)

# Output results
list(
  train_metrics = train_metrics,
  test_metrics = test_metrics
)
glm_test_accuracy <- test_metrics$accuracy

```




## 8. Model Evaluation

We evaluate the performance of each model using accuracy, which is the proportion of correct predictions over the total number of predictions.


## Output results

```{r}
#| label: model-evaluation
list( glm_accuracy = glm_test_accuracy,
  lda_accuracy = lda_accuracy,
  qda_accuracy = qda_accuracy,
  knn_accuracy = knn_accuracy,
  dt_accuracy = dt_accuracy
)

```


## Conclusion

We have successfully implemented various machine learning classification models to predict stock movement direction. The models were evaluated using accuracy, and further improvements can be made by tuning hyperparameters, adding more features, or using ensemble methods. 

Each model has its strengths and weaknesses, and the results from this tutorial can be used to guide further research or model enhancement.

