## 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(tidyverse)  # For data manipulation and visualization
library(caret)      # For model evaluation and cross-validation
library(rpart)      # For Decision Trees
library(Metrics)    # For additional evaluation metrics
library(janitor)    # For cleaning column names
library(MASS)
library(class)    # For KNN
library(boot)
```



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)
``` 


### 1. **Bootstrap Logistic Regression**
Apply bootstrap resampling to estimate logistic regression coefficients and evaluate variability.

```{r}
## Data cleaning and feature engineering


#| 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))

```


#### Define a Bootstrap Function:
```{r}
boot_logistic <- function(data, indices) {
  sampled_data <- data[indices, ]
  model <- glm(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, 
               data = sampled_data, family = binomial)
  return(coef(model))
}
```

#### Apply Bootstrapping:
```{r}
set.seed(123)
logistic_boot_result <- boot(data = psx_clean, statistic = boot_logistic, R = 1000)
print(logistic_boot_result)
```

#### Visualize and Extract Confidence Intervals:
```{r}
plot(logistic_boot_result)
```


```{r}
boot.ci(boot.out = logistic_boot_result, type = c("norm", "basic", "perc"))

```

---

### 2. **Bootstrap for LDA**

#### Define a Bootstrap Function:
```{r}
boot_lda <- function(data, indices) {
  sampled_data <- data[indices, ]
  lda_model <- lda(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, 
                    data = sampled_data)
  predictions <- predict(lda_model, newdata = data)$class
  accuracy <- mean(predictions == data$direction)
  return(accuracy)
}
```

#### Apply Bootstrapping: 

```{r}
set.seed(123)
lda_boot_result <- boot(data = psx_clean, statistic = boot_lda, R = 500)
print(lda_boot_result)
```

---

### 3. **Bootstrap for QDA**

#### Define a Bootstrap Function:
```{r}
boot_qda <- function(data, indices) {
  sampled_data <- data[indices, ]
  qda_model <- qda(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, 
                    data = sampled_data)
  predictions <- predict(qda_model, newdata = data)$class
  accuracy <- mean(predictions == data$direction)
  return(accuracy)
}
```

#### Apply Bootstrapping:
```{r}
set.seed(123)
qda_boot_result <- boot(data = psx_clean, statistic = boot_qda, R = 500)
print(qda_boot_result)
```

---

### 4. **Bootstrap for KNN**

#### Define a Bootstrap Function:
```{r}
boot_knn <- function(data, indices) {
  train_data <- data[indices, ]
  test_data <- data[-indices, ]
  knn_model <- 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
  )
  accuracy <- mean(knn_model == test_data$direction)
  return(accuracy)
}
```

#### Apply Bootstrapping:
```{r}
set.seed(123)
knn_boot_result <- boot(data = psx_clean, statistic = boot_knn, R = 500)
print(knn_boot_result)
```

---

### 5. **Bootstrap for Decision Trees**

#### Define a Bootstrap Function:
```{r}
boot_dt <- function(data, indices) {
  sampled_data <- data[indices, ]
  tree_model <- rpart(direction ~ lag1 + lag2 + lag3 + lag4 + lag5 + vol, 
                      data = sampled_data, method = "class")
  predictions <- predict(tree_model, newdata = data, type = "class")
  accuracy <- mean(predictions == data$direction)
  return(accuracy)
}
```

#### Apply Bootstrapping:
```{r}
set.seed(123)
dt_boot_result <- boot(data = psx_clean, statistic = boot_dt, R = 500)
print(dt_boot_result)
```

---

### Summarize Results:
Aggregate accuracy or coefficients across models to compare:
```{r}
list(
  logistic = mean(logistic_boot_result$t),
  lda = mean(lda_boot_result$t),
  qda = mean(qda_boot_result$t),
  knn = mean(knn_boot_result$t),
  decision_tree = mean(dt_boot_result$t)
)
```