R Notebook

library(tidyverse)
library(readxl)
library(corrplot)
#install.packages("Hmisc")
library(Hmisc)


set.seed(123)

setwd("C:/Users/Emily/Desktop/Brit School/2021 Fall/dow jones")

# read in data
dow_jones = read.csv("dow_jones_index.data")

# split into test and train
### Use quarter 1 (Jan-Mar) data for training and quarter 2 (Apr-Jun) data for testing.
train_dj = dow_jones %>% filter(quarter == 1)

test_dj = dow_jones %>% filter(quarter == 2)

head(train_dj)

head(test_dj)

nrow(train_dj)

[1] 360

nrow(test_dj)

[1] 390

train_dj %>%
  mutate(
    test_date = as.Date(date, format = "%m/%d/%Y")
  ) %>%
  group_by(test_date) %>%
  dplyr::summarize(count = n()) 

`summarise()` ungrouping output (override with `.groups` argument)

# count of nulls by column
train_dj %>%
  summarise_all(funs(sum(is.na(.))))

`funs()` is deprecated as of dplyr 0.8.0.
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

# visualizing time by each stock
train_dj %>%
  ggplot(aes(x = as.Date(date, format = "%m/%d/%Y"), y = percent_change_next_weeks_price)) +
  geom_line() +
  geom_point() +
  facet_wrap(~stock) +
  labs(title = "Percent Change Next Weeks Price")

train_dj %>%
  group_by(date) %>%
  dplyr::summarize(avg_chg = mean(percent_change_next_weeks_price, na.rm = TRUE)) %>%
  ggplot(aes(x = as.Date(date,, format = "%m/%d/%Y"), y  =avg_chg)) +
  geom_line() +
  labs(title = "Average Daily Percent Change - Next Weeks Price")

`summarise()` ungrouping output (override with `.groups` argument)

# correlations with variables that researchers used
# percent_change_price,
# percent_change_volume_over_last_wk, days_to_next_dividend, and percent_return_next_dividend

train_dj_subset = 
train_dj %>% 
  dplyr::select(percent_change_price, percent_change_volume_over_last_wk, previous_weeks_volume, percent_change_next_weeks_price, days_to_next_dividend, percent_return_next_dividend) %>%
  filter(!is.na(.))

head(train_dj_subset)

corrplot(
  cor(train_dj_subset, use = "pairwise.complete.obs"),
  type = "lower",
  addCoef.col = "black"
)

Decision Tree Model

Fitting a decision tree model using the subset of data only containing the variables the researchers used

# install and load needed packages
#install.packages("rpart")
library(rpart)

#install.packages('Metrics')
library(Metrics)

#install.packages('randomForest')
library(randomForest)

#install.packages('e1071')
library(e1071)

decision_tree_fit = 
  rpart(
    percent_change_next_weeks_price ~ percent_change_price + percent_change_volume_over_last_wk + previous_weeks_volume + days_to_next_dividend + percent_return_next_dividend,
    method = "anova",
    data = train_dj_subset
  )

# plot the tree
plot(decision_tree_fit,
     uniform = TRUE)
text(decision_tree_fit, use.n = TRUE, cex = .6)

dt_pred = predict(decision_tree_fit, test_dj, method = "anova")
rmse(actual = test_dj$percent_change_next_weeks_price, predicted = dt_pred)

[1] 3.144843

SVR

svr_fit = svm(    percent_change_next_weeks_price ~ percent_change_price + percent_change_volume_over_last_wk + previous_weeks_volume + days_to_next_dividend + percent_return_next_dividend,
                  data = train_dj_subset)

svr_pred = predict(svr_fit, test_dj)

rmse(predicted = svr_pred, actual = test_dj$percent_change_next_weeks_price)

[1] 2.945372

Tune SVR

svr_tuned_fit =tune(svm,
                    percent_change_next_weeks_price ~ percent_change_price + percent_change_volume_over_last_wk + previous_weeks_volume + days_to_next_dividend + percent_return_next_dividend,
                    data=train_dj_subset,
                    ranges=list(epsilon=seq(0,1,0.1), cost=  2^(2:9)))

#Print optimum value of parameters
print(svr_tuned_fit)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 6.7274 

plot(svr_tuned_fit)

NA
NA
NA

best_svr = svr_tuned_fit$best.model


best_svr_pred = predict(best_svr, test_dj)

rmse(predicted = best_svr_pred, actual = test_dj$percent_change_next_weeks_price)

[1] 2.972027

Now going to try the grid search on a smaller set of values since the darkest regions tend to be under cost = 20

svr_tuned_fit2 =tune(svm,
                    percent_change_next_weeks_price ~ percent_change_price + percent_change_volume_over_last_wk + previous_weeks_volume + days_to_next_dividend + percent_return_next_dividend,
                    data=train_dj_subset,
                    ranges=list(epsilon=seq(0,1,0.1), cost=  2^(2:5)))

#Print optimum value of parameters
print(svr_tuned_fit)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 6.7274 

plot(svr_tuned_fit2)

fit tuned SVR and calculate RMSE

best_svr2 = svr_tuned_fit2$best.model


best_svr_pred2 = predict(best_svr2, test_dj)

rmse(predicted = best_svr_pred2, actual = test_dj$percent_change_next_weeks_price)

[1] 2.972027

Trying the gridsearch on a smaller set of values did not help the RMSE over the first tuned model.

Results

The best performing model when comparing the decision tree fit to the SVR fit is the tuned SVR.