Bernauer_Andrew_final.Rmd

---
title: "Final Report"
author: "Andrew Bernauer"
date: "April 30, 2019"
output: html_document
---



#Difficult Learning Problems in Economics: Predicting Ether Price  









#a) Introduction 

For the semester long project I will be attempting to predict the known response variable ethereum price in terms of the predictors: block size, network hash rate growth rate, transactions, day of the week, and month. The project is regression based and not a classification task. Therefore the project will fall under the umbrella of supervised machine learning.    



$$\text{ether price} \approx \beta_{o} + \beta_{1} \times \text{blocksize} + \beta_{2} \times \text{networkhashgrowthrate} + \beta_{3}  \times \text{transactions} + \beta_{4} \times \text{day} + \beta_{5} \times \text{month} $$

Ethereum is a decentralized, open source, block chain technology, featuring smart contracts. The crypto currency that fuels the Ethereum block chain is Ether. Block size refers to the size of the Ethereum Block chain. Transactions are the number of transactions approved by the ledger. Hash rate is the rate at which it takes a miner to complete a Cryptographic computation on the block chain. 



 
```{r data carpentry, cache=FALSE}


library(scales)
library(ggplot2)
library(lubridate)
library(purrr)
library(dplyr)
library(readr)
library(tibble)
library(errorist)
library(broom)
library(caret)
library(tidyr)
library(randomForest)
library(gbm)

#substitute in the path on your machine to files 
ethereum_transaction_history <- read_csv("C:\\Users\\andre\\Documents\\ECON_490_ML\\ML_report 1\\ethereum-historical-data\\EthereumTransactionHistory.csv")

ethereum_block_size <-read_csv("C:\\Users\\andre\\Documents\\ECON_490_ML\\ML_report 1\\ethereum-historical-data\\EthereumBlockSizeHistory.csv"
)

ethereum_network_hash_rate <- read_csv("C:\\Users\\andre\\Documents\\ECON_490_ML\\ML_report 1\\ethereum-historical-data\\EthereumNetworkHashRateGrowthRate.csv")

ether_price <- read_csv("C:\\Users\\andre\\Documents\\ECON_490_ML\\ML_report 1\\ethereum-historical-data\\EtherPriceHistory(USD).csv")

# creating variable for day  
day_of_week <- 
  lubridate::mdy(ether_price$`Date(UTC)`) %>%
  wday( ,label=TRUE)  %>%
  sort()

# creating variable for month
month_of_year <-
  lubridate::mdy(ether_price$`Date(UTC)`) %>%
  month( ,label=TRUE) %>%
  sort()

signedlog10 <- function(x) {
  ifelse(abs(x) <= 1, 0, sign(x)*log10(abs(x)))
}

log_price <- signedlog10(ether_price[[3]])

log_history <- signedlog10(ethereum_transaction_history[[3]])

log_block_size <- signedlog10(ethereum_block_size[[3]])

log_hash_rate <- signedlog10(ethereum_network_hash_rate[[3]])

summary(log_price)

summary(log_history)

# construct list of variables to coerce into a tibble
l_ether <- list( price = ether_price[[3]],
                 transaction_history = as.integer(ethereum_transaction_history[[3]]),
                 block_size = as.integer(ethereum_block_size[[3]]),
                 hash_rate = ethereum_network_hash_rate[[3]], 
                 day = day_of_week, 
                 month = month_of_year,
                 date_utc = mdy(ethereum_network_hash_rate$`Date(UTC)`),
                 unix_time_stamp = ethereum_network_hash_rate$UnixTimeStamp,
                 log_10_price = log_price,
                 log_10_trans_history = log_history,
                 log_10_block_size = log_block_size,
                 log_10_hash_rate = log_hash_rate)

ether_df <- as_tibble(l_ether)


xs <-function(x){(x-mean(x))/(2*sd(x))}

```




```{r price_vs_time, warning=FALSE, cache=TRUE}
 
price_history <- ggplot(ether_df) +
  aes(date_utc, price) +
  geom_line(color="purple") + 
  xlab("Year")+
  ggtitle("Ethereum Price USD over time") + 
  scale_y_log10() 
  
  
  

 price_history
```


The price doesn't seem move at all significantly until 2016 so I decided to use a log transformation to the y axis. This removed a large amount the visual noise and made the price change overtime more evident. Besides that ethereum price starts to rise from around a dollar price to hovering around the ten dollars between 2016 and 2017. Peaks at a value of 1000 dollars in 2018.   





```{r block_size, cache=TRUE}
block_history <- ggplot(ether_df) +
  aes(date_utc, block_size) +
  geom_line(color ="purple") +
  scale_y_log10() +
  xlab("Year") +
  ylab("Block Size") +
  ggtitle("Ethereum Block Size Over Time")
  

block_history
```

Block size rises exponentially begining in 2017. Prior to this it fluctuates up and down between 1000 and 3000 which might be pointing to autocorrelation. In this same time period it peaks violently at 10000. Reaching a maximum of over 30000 in 2018. Another log tranformation was used in this plot.  



```{r network hashrate, cache=TRUE}
 hashrate_history <- ggplot(ether_df) + 
  aes(date_utc, hash_rate) +
  geom_line(color="purple") +
  scale_y_log10() +
  xlab("Year") +
  ylab("Network Hash rate Growth rate") +
  ggtitle("Ethereum Network Hash rate growth over time")

 hashrate_history
```

Network hashrate growth blows up exponentially and seems to follow that trend overtime. 

```{r transaction history, warning=FALSE, cache=TRUE}
 transaction_history <- ggplot(ether_df) + 
  aes(date_utc, transaction_history) +
  geom_line(color="purple") +
  xlab("Year") +
  ylab("Transactions") +
  scale_y_log10() +
  ggtitle("Ethereum Transactions over time")

 transaction_history
```

Transactions on the the Ethereum network are growing significantly overtime and more people are embracing it's use.





```{r bonus plot from eda, warning= FALSE, cache=TRUE}
 price_plot <- ggplot(ether_df) +
   aes(transaction_history, price) +
   geom_point(colour="purple") + 
   facet_wrap(~month) +
   geom_smooth(colour="yellow", alpha = 0.25) 

 price_plot +
   scale_x_log10("Transactions") +
   scale_y_log10("Ethereum Price USD($)") +
   ggtitle("Ethereum Price USD($) vs Transactions facetted by month") 
 
```

There appears to be some clustering of Ether price as transactions increases on the block chain increase isolating by month.



```{r}
library(ggplot2)
p <- ggplot(ether_df) +
  aes(x=log_10_hash_rate) +
  geom_histogram(bins = 34, fill = "purple") +
  ggtitle(label = "Histogram of the distribution of hash rate log 10 scale") +
  xlab("Distribution of log 10 transactions") + 
  ylab("Count")
  
  
p
```

Many of the $X_{1}, X_{2}, \dots, X_{n}$ or predictor variables in the dataset follow skewed distribution and feature multiple peaks. This could be pointing to the presence of sub-populations within the underlying data. 





```{r hist block_size,}
ggplot(ether_df) +
  aes(log_block_size) + geom_histogram(bins = 34, fill = 'purple') +
  xlab("Distribution of Ethereum BlockSize log 10 scale") +
  ylab("Count") +
  ggtitle("Histogram of Ethereum block size")

```

Likewise, blocksize has this similar property this sort of distribution doesn't suit itself well, for standard machine learning algorithms given the multiple peaks. Most of the variables follow gamma distributions in the cartesian coordinate scale. I declined to illustrate those graphs.


#b.) Summary Statistics 

```{r tidy_summary_code}

#summary stats for ethereum price  
ether_df %>% summarise(mean_price = mean(price, na.rm = TRUE), median_price = median(price), sd_price = sd(price), iqr_price = IQR(price), n = n(), mad_price = mad(price), min_price = min(price), max_price = max(price))

#summary stats for ethereum transaction history
ether_df %>% summarise(mean_transaction_history = mean(transaction_history), median_transaction_history = median(transaction_history), sd_transaction_history = sd(transaction_history), iqr_transaction_history = IQR(transaction_history), n = n(), mad_transaction_history = mad(transaction_history), min_transaction_history = min(transaction_history), max_transaction_history = max(transaction_history))

#summary stats for ethereum block size
ether_df %>% summarise(mean_block_size = mean(block_size), median_block_size = median(block_size), sd_block_size = sd(block_size), iqr_block_size = IQR(block_size), n = n(), mad_block_size = mad(block_size), min_block_size = min(block_size), max_block_size = max(block_size))

#summary stats for hash rate
ether_df %>% summarise(mean_hash_rate = mean(hash_rate), median_hash_rate = median(hash_rate), sd_hash_rate= sd(hash_rate), iqr_hash_rate = IQR(hash_rate), n = n(), mad_hash_rate = mad(hash_rate), min_hash_rate = min(hash_rate), max_hash_rate = max(hash_rate))




```








#c) Linear Regressions


```{r regressions,}
#running five regressions 
 
poly_reg <- lm(price ~ poly(transaction_history, 5) + poly(block_size, 5) + poly(hash_rate, 5), ether_df )


reg_obj <- lm(log(price) ~ log(block_size) + log(transaction_history) + log(hash_rate), data = ether_df,  subset = price > 0)

reg_obj_2 <- lm(log(price) ~ log(transaction_history) + log(hash_rate), data = ether_df, subset = price > 0)


reg_obj_3 <- lm(price ~ I(block_size)^2 + sqrt(transaction_history) + hash_rate, data = ether_df)

mols <- lm(price ~ block_size + transaction_history + hash_rate, data = ether_df)

log_ols <- lm(log(price) ~ log(transaction_history), data = ether_df, subset = price > 0)


 

```





The following code block returns the adjusted  $\hat{R}^{2}_{\text{adjusted}}$ statistics for the five regressions I ran.


```{r, cache=TRUE}
tidy_reg_obj <- reg_obj %>% 
   tidy()

tidy_reg_obj
```



```{r}
glance_tidy_poly <- glance(poly_reg)

glance_reg_obj <- glance(reg_obj)

glance_reg_obj_2 <- glance(reg_obj_2)

glance_reg_obj_3 <- glance(reg_obj_3)

glance_mols <- glance(mols)

glance_log_ols <- glance(log_ols)


glance_tidy_poly$adj.r.squared

glance_reg_obj$adj.r.squared

glance_reg_obj_2$adj.r.squared

glance_reg_obj_3$adj.r.squared

glance_mols$adj.r.squared

glance_log_ols$adj.r.squared

```

All of the coefficients return p values significant at the .10 cut off level while log block size is close. The $F$ statistic is significant; however the BIC and AIC below could be lower.










The log-log model with all the continuous predictors included had the highest adjusted $\,$ 
$\hat{R}^{2}_{\text{adjusted}} = 0.9363884$ $\,$ of all the regressions ran. A model excluding the log block size term comes in a close second. Another, option is the last model with just log transaction as a predictor. 

Predictors performing the best presented the following equation $$ \hat{y} = -10.65 + 0.072 \log(\text{Blocksize}) + 1.0007 \log(\text{TransactionHistory}) + 0.222\log(\text{Hashrate}) $$.

Which can be interpreted as a one percent increase in Blocksize resulting in a 7.2 percent increase in average Ether price holding all other variables equal, or fixed. The same interpretation applies to the other predictors with there respective coefficients.











```{r augment, cache=TRUE}
augmented_reg_obj <- reg_obj %>%
             augment()

some_plot <- ggplot(data = augmented_reg_obj, aes(log.price., .hat)) +
  geom_point(color = "purple", size = 1) +
  labs(x="log(price)", y="y hat") +
  ggtitle("Log price vs Y hat") 
  
  

some_plot
```



The preceding plot points to the non-linear nature of the data set. It may improve modeling results to include a polynomial term in future models.





```{r, warning=FALSE}
tidy_conf_it <- tidy(reg_obj, conf.int = TRUE)

confint_plot <- ggplot(tidy_conf_it, aes(term, estimate, color=term)) +
  geom_point() +
  geom_errorbar(aes(ymin=estimate - 1.96*std.error, ymax=estimate + 1.96*std.error)) 
  

confint_plot +
  ggtitle("Coefficients with standard errors") +
  scale_x_discrete()
  

```

The previous graph is misleading due to the scale with a scale adjustment it is more insightful.






```{r proper scale, warning=FALSE}

confint_plot + 
  scale_y_log10() +
  ggtitle("Coefficients with standard errors log-log scale") + 
  scale_x_discrete()

```





























```{r}
some_plot_2 <- ggplot(data = augmented_reg_obj, aes(.resid)) + 
  geom_histogram(binwidth = .3, fill = "purple") + 
  labs(x="Residuals", y="Frequency") +
  ggtitle("Distribution of Residuals")
  
  

some_plot_2
```
  




The majority of the residuals are distributed between -1 and 1 peaking at zero. This distribution is similar to the laplace distribution.  



```{r, warning= FALSE}
some_plot_3 <- ggplot(augmented_reg_obj, aes(log.price., .resid)) +
  geom_point(color = "purple", size = 1) +
  geom_hline(color = "yellow", yintercept = 0) 

some_plot_3 + 
  labs(x="Log Ether Price", y="Residuals")




```



The Residuals plotted against Ether Price suggest there is possible non-normality in the distribution of the errors.



















## d. and e.) Ridge and Lasso Regressions

```{r 1) 2) d) }
library(glmnet)
library(caret)


set.seed(385)


train_df <- ether_df %>% sample_frac(0.75)
test_df <- ether_df %>% setdiff(train_df)

train_df <- train_df[, -9, drop=FALSE]

test_df <- test_df[, -9, drop=FALSE]

covariates_train <- model.matrix(price~. , data = train_df)[, c(-1, -9), drop=FALSE]

covariates_test <- model.matrix(price~. , data=test_df)[, c(-1, -9), drop=FALSE]
response_train <- as.numeric(unlist(train_df[[1]]))
response_test <- as.numeric((unlist(test_df[[1]])))

grid <- 10^seq(5, -5, length.out = 100)

ridge_model <- glmnet(covariates_train, response_train, alpha = 0, family ="gaussian")

plot.glmnet(ridge_model, xvar = "lambda")

lasso_model <- glmnet(covariates_train, response_train, alpha = 1,  family = "gaussian")


plot.glmnet(lasso_model, xvar = "lambda", label = "true")




```

A $\lambda = 49.75$, for the hyperparameter of the ridge model significantly penalizes the coefficients of the model. Shrinking all of the coefficients close to zero displayed in the plot above. 

Additionally, the lasso model has a $\lambda = 2.26$ which ends up culling eleven of the twenty five coefficients assocaited with the model matrix. Which is evident in the plot above. 




```{r 2a)}
set.seed(385)

cv_lasso <- cv.glmnet(covariates_train, train_df$price, alpha=1)


plot(cv_lasso)

set.seed(385)
cv_ridge <- cv.glmnet(covariates_train, train_df$price, alpha=0)

plot(cv_ridge)

best_lambda_ridge <- cv_ridge$lambda.1se

best_lambda_lasso <- cv_lasso$lambda.1se

print(best_lambda_ridge)

print(best_lambda_lasso)




```

I elect to pick a lambda within one standard deviation as a way to avoid overfitting the model to the training data and performing poorly on the test data. 








```{r 2) e)}
ridge_coefficients <- coef.glmnet(cv_ridge, newx = covariates_train,  s = best_lambda_ridge)



lasso_coef <- coef.glmnet(cv_lasso, newx = covariates_train, s = best_lambda_lasso)

print(lasso_coef)

print(ridge_coefficients)
```

Lasso as an operator ending up culling the following coefficients: day.L, day.Q, day.C, day^4, Month.L, date_utc, unix_timestamp, and the other coefficients for transaction history, block size, hash rate in the log 10 scale. The model coefficients are not the most interpratable as they were produced via the model matrix 

The same can be said for the ridge regression model.





```{r 2) g),}
lasso_test_pred = predict.glmnet(lasso_model, newx = covariates_test, s=best_lambda_lasso)

ridge_test_pred = predict.glmnet(ridge_model, newx = covariates_test, s=best_lambda_ridge)

lasso_mse <- mean((response_test - lasso_test_pred) ^ 2)

ridge_mse <- mean((response_test - ridge_test_pred) ^ 2)
```

The **Lasso ** model returned an $MSE$ of 88166.71 while the **Ridge** model returned a more modest $MSE$ of 6682.426 on the test sets respectively. The Ridge regression model returns the lowest $MSE $of the two methods though the performance is not great.
 




##f.) Regression trees 

```{r 3) a-d, }
library(tree)
set.seed(385)

ether_tree <- tree(price~., data = train_df)

cv_tree <- cv.tree(ether_tree, K=10 )
plot(cv_tree$size, cv_tree$dev, type = 'b')

pruned_tree_df <- prune.tree(ether_tree, best = 5)
plot(pruned_tree_df)
text(pruned_tree_df, pretty = 0)

treee_fit = predict(pruned_tree_df, newdata = test_df)


single_tree_fit <- mean((response_test - treee_fit) ^ 2)

single_tree_fit



```

Variables used in the final pruned tree are: transaction history which represent the number of transactions on the Ethereum block chain and unix time stamp which is a running count of the seconds since Jan 01 1970. The root of the tree includes a split with the variable transaction history splitting on the condition that if transaction history is less than one hundred sixty two thousand six hundred and eighty three transactions. After that split follows, another if then style statement with $\text{transaction history} < 162,683$ yielding an average ether price prediction of 13.83. On the other hand if that condition is not satisfied average ether price yields a prediction of 288.50. The splits on the other side of the tree include a split, for $\text{transaction history} < 876,929 $, and if $\text{transaction history} > 867,929$ yields an average ether price prediction of $ 960.40$. Under the split above is another split under unix time which essentially prefers ether price fall over time.  

The pruned regression tree contains a total of five terminal nodes each displaying the average prediction based on those splits in the tree. The splits function to partition space into mutually exclusive regions. 











The $MSE$ of the bootstrapped trees is a performance upgrade over the results from the previous performance of the single regression tree trained on the training data. This is not completely surprising as we would expect the model to perform better on our training data and to underperform on the test data. Surprisingly though I was able to get an  $MSE$ that was significanctly better accross multiple seeds and fitting trees which is promising, for future reports.  Of course there is a trade-off in that we don't want our model to overfit the training data and not perform well on the test set. 












```{r #1, #2}
set.seed(385)

new_train_ether <- ether_df[,  -9, drop=FALSE] %>%
  sample_frac(0.7) 
  
new_test_ether <- ether_df[,  -9, drop=FALSE] %>%
  setdiff(new_train_ether)

bagged.ether <- randomForest(price ~ ., data = new_train_ether,
                             ntree=400, mtry=4, importance=TRUE)
bagged.ether



```


# g.) Bagging


```{r}
plot(bagged.ether)
```


The Errors start to level off as trees start approaching the 50 to 75 tree range. Increasing the number of trees does not impact the error rate that much, as the number of trees gets large. 






```{r}
yhat_bagged <- predict(bagged.ether, newdata = new_test_ether)

ggplot() + geom_point(aes(x = yhat_bagged, y = price), data = new_test_ether, color = "purple", alpha = 0.50) + geom_abline(color = "yellow") + 
  labs(x="Bagging Estimation of Ether", y="Ether test Price", title="Bagging Estimation of Ether vs Actual Test Ether Price") 
  
```

The Bagging estimation of Ether Price presents some of the best predictions on the test set of all methods. However the fit is not as tight as Ether Price increases. In addition it is an improvement over the bootstrapped version of the trees from my previous report. 




```{r}
bagging_mse <- round(mean((new_test_ether$price - yhat_bagged)^2), 2)
bagging_mse
lasso_mse
ridge_mse
```
The bagging model is the best $MSE = 1023.55$; compared to the large $MSE =  5159.674$ for lasso and the large mse of $MSE = 6704.901$  for ridge regression. Meaning that the model has less variance in performance about the mean. 





```{r, 2e)}
importance(bagged.ether)
```




```{r}
varImpPlot(bagged.ether)
```

In terms of important variables for the bagging model date, unix time stamp, log_10_transaction history, and transaction history have the largest increased percentage in the mean squared errors. It is a more robust measure of error versus the IncNodePurity. 


# h.) Random Forest
```{r}
set.seed(385)

ether_forest <- randomForest(price ~ ., data = new_train_ether,
                             mtry = 10,
                             importance = TRUE,
                             do.trace = 100)
plot(ether_forest)
```

Error mostly flat lines as number of trees approaches 50.


```{r}
y_hat_rf <- predict(ether_forest, newdata = new_test_ether)

```



```{r}
# Set mtry using hyperparamter tuning

oob_err <- double(10)
test_err <- double(10)

#mtry is no of Variables randomly chosen at each split
for(mtry in 1:10) 
{
  rf_ether <- randomForest(price ~ . , data = new_train_ether, mtry=mtry, ntree=400) 
  oob_err[mtry] <- rf_ether$mse[400] #Error of all Trees fitted on training
  
  pred <- predict(rf_ether, new_test_ether) #Predictions on Test Set for each Tree
  test_err[mtry] <- with(new_test_ether, mean( (price - pred)^2)) # "Test" Mean Squared Error
  
  print(mtry)
  
}

```

```{r}
round(test_err, 2)
```
```{r}
round(oob_err, 2)
```
An mtry of four returns the lowest out of bag error for random forest. 



```{r}
rf.opt <- randomForest(price ~ ., data = new_train_ether, n_tree = 400, mtry = 4)

y_hat_rf_opt <- predict(rf.opt, newdata = new_test_ether)

ggplot() + geom_point(aes(x = y_hat_rf_opt, y = new_test_ether$price), color = "purple", alpha=0.5) + geom_abline(color = "yellow") + 
  labs(x="Random Forest estimation of Ether Price", y="Ether Price", title = "Random Forest estimation of Ether Price vs Actual Ether Price")

```


The residuals increase as Ether Price increases similar to other tree based methods, but has one of the best fits of all the tree based models.




```{r}
mse_random_forest <- round(mean((new_test_ether$price - y_hat_rf_opt)^2))
lasso_mse
ridge_mse
bagging_mse
mse_random_forest
```
Random forest fits Ether Price really well better than any other models used so far, except bagging. 



```{r}
varImpPlot(rf.opt)
```
Like the other tree based models random forest emphasizes the importance of transaction history as well as time.



```{r}
matplot(1:mtry , cbind(oob_err,test_err), pch=20 , col=c("red","blue"),type="b",ylab="Mean Squared Error",xlab="Number of Predictors Considered at each Split")
legend("topright",legend=c("Out of Bag Error","Test Error"),pch=19, col=c("red","blue"))
```
   




The ideal number of predictors Considered at each split is different depending on whether our goal is to minimize the mean squared error of the test error or the out of bag error the best number of predictors to consider at each split happens to be 4 predictors. While 7 predictors allow for a minimized mean squared error. However; for cross validation we will use the ideal number of predictors based on the out of bag error. 

# i.) Boosting

```{r}
set.seed(385)
boosting_ether <- gbm(price ~ . -date_utc, 
                      data = new_train_ether,
                      distribution = "gaussian",
                      n.trees = 5000,
                      interaction.depth = 4) 
summary(boosting_ether)
```

Transaction history, Unix time stamp, and hash rate are all displayed as important variables according to the importance matrix. Transaction history carries the largest relative importance of all the variables at a value of 76.36 pecent. While unix time stamp has an importance of 11.67 percent which is not very significant. The only other variable containing a significant amount of importance is network hash rate. 



```{r}
yhat_boost <- predict(boosting_ether, newdata = new_test_ether, n.trees = 5000)

ggplot() + geom_point(aes(x = yhat_boost, y = new_test_ether$price), color = "purple", alpha=0.5) + geom_abline(color = "yellow") + 
  labs(x="Boosting Estimation of Ether Price", y="Ether Price test set", title = "Boosting Estimation of Ether Price vs Actual Ether Price")
```

Boosting is a relatively viable modeling technique, for predicting ether price given the fit above.







```{r}
boosting_mse <- mean((new_test_ether$price - yhat_boost)^2)
boosting_mse
lasso_mse
ridge_mse
bagging_mse
```

Boosting returns a very low Mean Square error of 1023.55.

## j.) XGboost

```{r}
library(xgboost)
f <- price ~ transaction_history + block_size + hash_rate + day + month + unix_time_stamp + log_10_trans_history + log_10_block_size + log_10_hash_rate

X_training <- model.matrix(f, data = new_train_ether)

Y_training <- as.matrix(new_train_ether$price)

d_train <- xgb.DMatrix(data = X_training, label = Y_training)

X_test <- model.matrix(f, data = new_test_ether)

set.seed(385)
ether.xgb = xgboost(data=d_train,
                     max_depth=5,
                     eta = 0.1,
                     nrounds=40, 
                     lambda=0,
                     print_every_n = 10,
                     objective="reg:linear")



```



```{r}
yhat_xgboost <- predict(ether.xgb, newdata = X_test)

ggplot() + geom_point(aes(x = yhat_xgboost, y = price), data = new_test_ether, color = "purple", alpha = 0.25) + 
  geom_abline(color = "yellow") + 
  labs(x="Xgboost estimation of ether price", y="Actual ether price", title = "Xgboost estimation of Ether price vs Actual Test ether price") 
  
```
The xgboosted model performs relatively strong, but like other tree based models the fit is not as good as Ether Price increases. 




```{r}
 xg_boosted_mse <- mean((new_test_ether$price - yhat_xgboost)^2) 
 lasso_mse
 ridge_mse
 bagging_mse
 mse_random_forest
 boosting_mse
 xg_boosted_mse
```

An MSE of 1415.984 is returned by the xgboosted model is respectable in comparison to boosting and the other tree-based methods, but does not perform aswell as Random Forest or gradient boosting. 







```{r}
 importance_xgb <- xgb.importance(colnames(X_training), model = ether.xgb)

 importance_xgb
 
 xgb.plot.importance(importance_matrix = importance_xgb, rel_to_first = TRUE, xlab="Relative Importance")
```

Transaction history displays an overwhelming level of relative importance for the xgboosted model as shown above. With time or unix time stamp being the only other important variable.


```{r}
library(caret)

f <- price ~ transaction_history + block_size + hash_rate + day + month + unix_time_stamp + log_10_trans_history + log_10_block_size + log_10_hash_rate

X_training <- model.matrix(f, data = new_train_ether)

Y_training <- as.matrix(new_train_ether$price)

d_train <- xgb.DMatrix(data = X_training, label = Y_training)

X_test <- model.matrix(f, data = new_test_ether)

params_xgboost <- list(nrounds = seq(40, 250, 5),
                     max_depth = 4,
                     eta = seq(2, 30,5)/60,
                     gamma = 1,
                     colsample_bytree = 1,
                     min_child_weight = 9,
                     subsample =  1.0
                     )

hyper_grid_xgboost <- expand.grid(params_xgboost)
xgb_control <- trainControl(method="cv",
                            number = 3)
set.seed(385)
xgboost_tuned <- train(price~.,
                        data=new_train_ether,
                        trControl=xgb_control,
                        tuneGrid=hyper_grid_xgboost,
                        lambda=0,
                        method="xgbTree")
xgboost_tuned$bestTune
```

An $\eta = 0.2$ implies the model starts out as a slow learner and requires a larger number of boosting iterations to perform well, but is more robust to errors as a result.  


```{r}
plot(xgboost_tuned)
```
The previous plot shows this as the case. 

```{r}
tuned.xgb <- xgb.train(params = xgboost_tuned$bestTune,
                       data = d_train, 
                       nrounds = 225)
yhat_xgboost_tuned <- predict(tuned.xgb, X_test)
tuned_xgboost_mse <- mean((yhat_xgboost_tuned - new_test_ether$price)^2)
tuned_xgboost_mse
xg_boosted_mse
```
The hyper-paramater tuned model only offers a marginal decrease in the $MSE$ aproximately 56. 


```{r}
gbm.grid <- expand.grid(interaction.depth = c(1, 4, 5, 9),
                        n.trees = (1:100)*50,
                        n.minobsinnode = c(10, 20, 25),
                        shrinkage = 0.1)
nrow(gbm.grid)

fit.control <- trainControl(method = "repeatedcv",
                            number = 5,
                            repeats = 5)

set.seed(385)

gbm.tuned <- train(price ~ . -date_utc, 
                   data = new_train_ether,
                   method = "gbm",
                   trControl = fit.control,
                   verbose = FALSE,
                   tuneGrid = gbm.grid)

gbm.tuned




```



```{r}
y_hat_tuned <- predict.train(gbm.tuned, newdata = new_test_ether, type ='raw')

ggplot()+ geom_point(aes(x = y_hat_tuned, y = new_test_ether$price), color = "purple", alpha = 0.50) + 
  geom_abline(color = "yellow") + 
  labs(x="Parameter tuned boosting estimate of Ether Price", y="Actual test-set Price", title = "Parameter tuned boosting estimation of Ether Price vs Actual test-set Price")
```

The Parameter boosted model gives a relatively tight fit around the ideal fit line. With larger residuals resulting as the price increases. 


```{r}
mse_tuned <- round(mean((new_test_ether$price - y_hat_tuned)^2), 2)
mse_tuned
mse_random_forest
bagging_mse
lasso_mse
ridge_mse
boosting_mse
xg_boosted_mse
```
The parameter tuned boosting model returns the third lowest mean squared errors of  any of the models. Though I should be cautious of overfitting as I grew a large number of trees. 


```{r, cache=TRUE}
with_2_standard_errs <- oneSE(gbm.tuned$results, metric = "RMSE", num = 5, maximize = FALSE)
gbm.tuned$results[with_2_standard_errs, 1:6]
```
These would be the ideal hyperparameter values to fit a gradient boosted model to avoid overfitting based on the one standard error rule of thumb. 


```{r, cache=TRUE}
summary(gbm.tuned)
```

Similiar to previous importance plots transaction history has the largest relative influence on ether price. Which makes some intuitive sense as investors might believe in the technology more as it usage increases. None of the other variables show much of a relative influence. 


Random forest appears to be the second best model, for modeling Ether price it avoids overfitting. In addition, the model constructed used the out of bag error to determine the number of splits to consider, four predictors. So we could expect the the results to generalize well to new test data. Which is not the case for the parameter tuned gbm model having higher mean squared error. 

#k) Neural Network
```{r}
library(keras)
set.seed(385)
 keras_train <- ether_df %>%
   sample_frac(.75)
 keras_test <- ether_df %>%
   setdiff(keras_train)
 keras_train_data <- as.matrix(keras_train[!names(ether_df) %in% c("price","log_10_price", "day", "month", "date_utc")])
 keras_test_data <- as.matrix(keras_test[!names(ether_df) %in% c("price","log_10_price", "day", "month", "date_utc")])
 keras_train_labels <- as.matrix(keras_train[, "price"])
 keras_test_labels <- as.matrix(keras_test[, "price"])
 
 keras_train_data <- scale(keras_train_data)
 keras_test_data <- scale(keras_test_data)


 
 
build_model <- function(){
 model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation="relu",
                input_shape = dim(keras_test_data)[2]) %>%
     layer_dense(units = 64, activation = "sigmoid") %>%
     layer_dense(units = 64, activation = "exponential") %>%
     layer_activation_leaky_relu(alpha = 0.01) %>%
     layer_dense(units = 1) 
 
 model %>% compile(
      loss = "mse",
    optimizer = optimizer_rmsprop(),
      metrics = list("mean_absolute_error"))
 
 
model
}
model <- build_model()
model %>% summary()
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 68)
epochs <- 325
# set the X and Y 

history <- model %>% fit(
  keras_train_data,
  keras_train_labels,
  epochs = epochs,
  validation_split = 0.50
)

```


```{r}
library(ggplot2)
plot(history, metrics = "mean_absolute_error", smooth = FALSE) 
```


The mean absolute error stabilizes around 100 epochs or iterations 



```{r}
model_2 <- build_model()

model_2 %>% compile(
  loss="mse",
  optimizer = optimizer_rmsprop(lr = 0.022),
  metrics = list("mean_absolute_error")
)
history2 <- model_2 %>% fit(
  keras_train_data,
  keras_train_labels,
  epochs=100,
  validation_split = 0.50,
  callbacks = list(early_stop)
)
  
```

```{r}
plot(history2, metrics = 'mean_absolute_error', smooth = FALSE)
```

The mean absolute errors are not as pronounced in the tuned model.

```{r}
nueral_predictions <- model %>% predict(keras_test_data)
net_mse <- mean((keras_test_labels - nueral_predictions)^2)
net_mse
```

A Neural Network performs decently given the complexity of the dataset. I experimented with different activation functions and numbers of layers. Finding success with interweaving layers of relu activation functions as well as exponential activiation functions. This greatly increased the performance of the model though I started having issues where the errors were not converging and took that as a sign that the vanishing gradient problem was taking place. One solution is to include a leaky relu function because it has a negative slope as x decreases it will still register a neuron as firing while a relu won't. This increased the performance of my model more significantly than adding more layers.


#i.) extra credit


```{r knn regression, eval=FALSE, cache=TRUE}
library(caret)

set.seed(20)

train_index <- createDataPartition(ether_df$price, p = .8, list = FALSE, times = 1)

head(train_index)

ether_df_train <- ether_df[train_index, , drop = FALSE]

ether_df_test <- ether_df[-train_index, , drop = FALSE]

tr_control <- trainControl(method = "repeatedcv",
             number = 10,
             repeats = 3
             )
knn_model <-  train(price ~., data = ether_df,
                    method = "knn",
                    preProcess = c("center", "scale"), 
                    trControl = tr_control, 
                    metric = "RMSE",
                    tuneLength = 10
                    )

knn_model

knn_model$results

y_hat <- predict(knn_model, newdata = ether_df)

predicted <- predict(knn_model, ether_df)

plot(ether_df$price, predicted)
sqrt(sum((predicted - ether_df$price) ^ 2) / length(ether_df)) 

```

```{r 4)a-c, cache=TRUE}
set.seed(385)

statistic <- function(x, n=100, best=5){
  tree_1 <- tree(price~., data = train_df)
  pruned_tree_1 <- prune.tree(tree_1, k=10, best = best)
  
 return(  predict(pruned_tree_1, newdata = test_df) )
  
  
  
}

 result <- boot::boot(data = train_df, statistic, R=100, sim = "ordinary")
 
mean((test_df$price - result$t0)^ 2)


```

#m.) comparing all models 

```{r}

single_tree_fit
lasso_mse
ridge_mse
bagging_mse
mse_random_forest
boosting_mse
mse_tuned
xg_boosted_mse
tuned_xgboost_mse
net_mse


```

The lowest $MSE$ was produced by the bagging model, followed by Random Forest, hyper-parameter tuned boosting, boosting, hyper-parameter tuned xgboost, xgboost, the neural network, lasso regression, ridge regression, and finally regression trees. 



| Machine learning method | MSE     |
|:-----------------------:|---------|
| Bagging                 | 1023.55 |
| Random Forest           | 1035    |
| Tuned boosting          | 1172.63 |
| Boosting                | 1299.22 |
| Tuned xgboost           | 1360.19 |
| Xgboost                 | 1415.98 |
| Neural Network          | 3596.69 |
| Lasso                   | 5159.67 |
| Ridge                   | 6704.90 |
| Regression Tree         | 8066.87 |


#n.) Conclusion 

Of all the machine learning techniques covered bagging proved the most accurate model for predicting Ether price. While tree methods generally performed better than other machine learning methods.  The techniques covered in class don't lend themselves naturally to dealing with time series data if we had covered Neural Networks in further depth I believe a neural network could perform better than all the other models. As neural networks lend themselves well with dealing with sequence data and in particular time series data one such technique is LSTM or Long short-term memory. Another possible approach might consider using lags which is out of the scope of the class, but is useful in time series modeling.