主要議題：預測股票的投資報酬
學習重點：
rm(list=ls(all=T))
Sys.setlocale("LC_ALL","C")## [1] "C"options(digits=5, scipen=12)
library(dplyr)
library(caTools)
library(caret)
library(flexclust)Load StocksCluster.csv into a data frame called “stocks”.
stocks = read.csv("data/StocksCluster.csv")How many observations are in the dataset? + 11580
What proportion of the observations have positive returns in December?
mean(stocks$PositiveDec)## [1] 0.54611What is the maximum correlation between any two return variables in the dataset? You should look at the pairwise correlations between ReturnJan, ReturnFeb, ReturnMar, ReturnApr, ReturnMay, ReturnJune, ReturnJuly, ReturnAug, ReturnSep, ReturnOct, and ReturnNov.
cor(stocks[,1:11]) %>% round(2) ## ReturnJan ReturnFeb ReturnMar ReturnApr ReturnMay ReturnJune
## ReturnJan 1.00 0.07 -0.09 -0.04 -0.04 0.09
## ReturnFeb 0.07 1.00 -0.16 -0.19 -0.10 0.17
## ReturnMar -0.09 -0.16 1.00 0.01 0.00 -0.09
## ReturnApr -0.04 -0.19 0.01 1.00 0.06 -0.01
## ReturnMay -0.04 -0.10 0.00 0.06 1.00 -0.02
## ReturnJune 0.09 0.17 -0.09 -0.01 -0.02 1.00
## ReturnJuly -0.08 -0.06 0.00 0.08 0.09 -0.03
## ReturnAug -0.02 0.13 -0.02 -0.05 -0.03 0.01
## ReturnSep -0.03 0.04 0.08 -0.03 0.02 0.04
## ReturnOct 0.14 -0.09 -0.01 0.05 0.02 -0.02
## ReturnNov 0.07 -0.15 0.04 0.03 0.05 -0.07
## ReturnJuly ReturnAug ReturnSep ReturnOct ReturnNov
## ReturnJan -0.08 -0.02 -0.03 0.14 0.07
## ReturnFeb -0.06 0.13 0.04 -0.09 -0.15
## ReturnMar 0.00 -0.02 0.08 -0.01 0.04
## ReturnApr 0.08 -0.05 -0.03 0.05 0.03
## ReturnMay 0.09 -0.03 0.02 0.02 0.05
## ReturnJune -0.03 0.01 0.04 -0.02 -0.07
## ReturnJuly 1.00 0.00 0.07 -0.05 -0.05
## ReturnAug 0.00 1.00 0.00 -0.08 -0.12
## ReturnSep 0.07 0.00 1.00 -0.06 -0.02
## ReturnOct -0.05 -0.08 -0.06 1.00 0.19
## ReturnNov -0.05 -0.12 -0.02 0.19 1.00colMeans(stocks[,1:11]) %>%sort # %>% barplot(las=2, cex.names=0.8, cex.axis=0.8)## ReturnSep ReturnFeb ReturnJuly ReturnOct ReturnJune ReturnNov
## -0.0147208 -0.0076048 0.0030509 0.0056508 0.0059379 0.0113874
## ReturnJan ReturnAug ReturnMar ReturnMay ReturnApr
## 0.0126316 0.0161983 0.0194023 0.0247366 0.0263081Which month (from January through November) has the largest mean return across all observations in the dataset?
Which month (from January through November) has the smallest mean return across all observations in the dataset?
Run the following commands to split the data into a training set and testing set, putting 70% of the data in the training set and 30% of the data in the testing set:
set.seed(144)
spl = sample.split(stocks$PositiveDec, SplitRatio = 0.7)
stocksTrain = subset(stocks, spl == TRUE)
stocksTest = subset(stocks, spl == FALSE)
set.seed(144)
spl = sample.split(stocks$PositiveDec, SplitRatio = 0.7)
stocksTrain = subset(stocks, spl)
stocksTest = subset(stocks, !spl)Then, use the stocksTrain data frame to train a logistic regression model (name it StocksModel) to predict PositiveDec using all the other variables as independent variables. Don’t forget to add the argument family=binomial to your glm command.
stocksmodel= glm(PositiveDec~., stocksTrain, family = "binomial")
pred= predict(stocksmodel, type="response" )
table(stocksTrain$PositiveDec, pred>=0.5)##
## FALSE TRUE
## 0 990 2689
## 1 787 3640(990+3640)/nrow(stocksTrain)## [1] 0.57118What is the overall accuracy on the training set, using a threshold of 0.5? + 0.567
pred2= predict(stocksmodel, newdata= stocksTest, type="response" )
table(stocksTest$PositiveDec, pred2>=0.5)##
## FALSE TRUE
## 0 417 1160
## 1 344 1553(417+1553)/nrow(stocksTest)## [1] 0.56707Now obtain test set predictions from StocksModel. What is the overall accuracy of the model on the test, again using a threshold of 0.5? + 0.54606
mean(stocksTest$PositiveDec) ## [1] 0.54606What is the accuracy on the test set of a baseline model that always predicts the most common outcome (PositiveDec = 1)?
Now, let’s cluster the stocks. The first step in this process is to remove the dependent variable using the following commands:
limitedTrain = stocksTrain
limitedTrain$PositiveDec = NULL
limitedTest = stocksTest
limitedTest$PositiveDec = NULL
#LTR= stocksTrain[,1:11]
#LTS= stocksTest[,1:11]Why do we need to remove the dependent variable in the clustering phase of the cluster-then-predict methodology? + Needing to know the dependent variable value to assign an observation to a cluster defeats the purpose of the methodology
In the market segmentation assignment in this week’s homework, you were introduced to the preProcess command from the caret package, which normalizes variables by subtracting by the mean and dividing by the standard deviation.
In cases where we have a training and testing set, we’ll want to normalize by the mean and standard deviation of the variables in the training set. We can do this by passing just the training set to the preProcess function:
preproc = preProcess(limitedTrain)
normTrain = predict(preproc, limitedTrain)
normTest = predict(preproc, limitedTest)
mean(normTrain$ReturnJan)## [1] 2.1006e-17mean(normTest$ReturnJan)## [1] -0.00041859What is the mean of the ReturnJan variable in normTrain?
What is the mean of the ReturnJan variable in normTest?
Why is the mean ReturnJan variable much closer to 0 in normTrain than in normTest?
Set the random seed to 144 (it is important to do this again, even though we did it earlier). Run k-means clustering with 3 clusters on normTrain, storing the result in an object called km.
set.seed(144)
km = kmeans(normTrain, 3)
table(km$cluster)##
## 1 2 3
## 3157 4696 253Which cluster has the largest number of observations?
Recall from the recitation that we can use the flexclust package to obtain training set and testing set cluster assignments for our observations (note that the call to as.kcca may take a while to complete):
library(flexclust)
km.kcca = as.kcca(km, normTrain)## Found more than one class "kcca" in cache; using the first, from namespace 'kernlab'## Also defined by 'flexclust'## Found more than one class "kcca" in cache; using the first, from namespace 'kernlab'## Also defined by 'flexclust'clusterTrain = predict(km.kcca)## Found more than one class "kcca" in cache; using the first, from namespace 'kernlab'
## Also defined by 'flexclust'clusterTest = predict(km.kcca, newdata=normTest)
table(clusterTest)## clusterTest
## 1 2 3
## 1298 2080 96How many test-set observations were assigned to Cluster 2?
Using the subset function, build data frames stocksTrain1, stocksTrain2, and stocksTrain3, containing the elements in the stocksTrain data frame assigned to clusters 1, 2, and 3, respectively (be careful to take subsets of stocksTrain, not of normTrain). Similarly build stocksTest1, stocksTest2, and stocksTest3 from the stocksTest data frame.
#stocksTrain1 = subset(stocksTrain, clusterTrain==1)
#stocksTrain2 = subset(stocksTrain, clusterTrain==2)
#stocksTrain3 = subset(stocksTrain, clusterTrain==3)
#stocksTest1 = subset(stocksTest, clusterTrain==1)
#stocksTest2 = subset(stocksTest, clusterTrain==2)
#stocksTest3 = subset(stocksTest, clusterTrain==3)
tapply(stocksTrain$PositiveDec , clusterTrain ,mean)## 1 2 3
## 0.60247 0.51405 0.43874tapply(stocksTest$PositiveDec , clusterTest ,mean)## 1 2 3
## 0.61402 0.51250 0.35417Which training set data frame has the highest average value of the dependent variable?
Build logistic regression models StocksModel1, StocksModel2, and StocksModel3, which predict PositiveDec using all the other variables as independent variables. StocksModel1 should be trained on stocksTrain1, StocksModel2 should be trained on stocksTrain2, and StocksModel3 should be trained on stocksTrain3.
#StocksModel1= glm(PositiveDec~., stocksTrain1, family = "binomial")
#StocksModel2= glm(PositiveDec~., stocksTrain2, family = "binomial")
#StocksModel3= glm(PositiveDec~., stocksTrain3, family = "binomial")
#summary(StocksModel1)
#summary(StocksModel2)
#summary(StocksModel3)
M = lapply(split(stocksTrain, clusterTrain), function(x)
glm(PositiveDec~., data=x, family="binomial") )
sapply(M, function(x) coef(summary(x))[,1])## 1 2 3
## (Intercept) 0.172240 0.10293 -0.1818958
## ReturnJan 0.024984 0.88451 -0.0097893
## ReturnFeb -0.372074 0.31762 -0.0468833
## ReturnMar 0.595550 -0.37978 0.6741795
## ReturnApr 1.190478 0.49291 1.2814662
## ReturnMay 0.304209 0.89655 0.7625116
## ReturnJune -0.011654 1.50088 0.3294339
## ReturnJuly 0.197692 0.78315 0.7741644
## ReturnAug 0.512729 -0.24486 0.9826054
## ReturnSep 0.588327 0.73685 0.3638068
## ReturnOct -1.022535 -0.27756 0.7822421
## ReturnNov -0.748472 -0.78747 -0.8737521Which variables have a positive sign for the coefficient in at least one of StocksModel1, StocksModel2, and StocksModel3 and a negative sign for the coefficient in at least one of StocksModel1, StocksModel2, and StocksModel3? Select all that apply.
Using StocksModel1, make test-set predictions called PredictTest1 on the data frame stocksTest1. Using StocksModel2, make test-set predictions called PredictTest2 on the data frame stocksTest2. Using StocksModel3, make test-set predictions called PredictTest3 on the data frame stocksTest3.
#PredictTest1= predict(StocksModel1, newdata=stocksTest1, type="response")
#PredictTest2= predict(StocksModel1, newdata=stocksTest2, type="response")
#PredictTest3= predict(StocksModel1, newdata=stocksTest3, type="response")
Pred = lapply(1:3, function(x)
predict(M[[x]], stocksTest[clusterTest==x,], type='response') ) #告訴它真的是x，不是第欄。而X是指定的1:3
sapply(1:3, function(x)
table(stocksTest$PositiveDec[clusterTest==x], Pred[[x]] > 0.5) %>% {sum(diag(.))/sum(.)} )## [1] 0.61941 0.55048 0.64583What is the overall accuracy of StocksModel1 on the test set stocksTest1, using a threshold of 0.5?
#table(stocksTest1$PositiveDec, PredictTest1>0.5)What is the overall accuracy of StocksModel2 on the test set stocksTest2, using a threshold of 0.5?
#table(stocksTest2$PositiveDec, PredictTest2>=0.5)What is the overall accuracy of StocksModel3 on the test set stocksTest3, using a threshold of 0.5?
#table(stocksTest3$PositiveDec, PredictTest3>=0.5)To compute the overall test-set accuracy of the cluster-then-predict approach, we can combine all the test-set predictions into a single vector and all the true outcomes into a single vector:
table( do.call(c, split(stocksTest$PositiveDec,clusterTest)), do.call(c, Pred) > 0.5 ) %>%
{sum(diag(.))/sum(.)}## [1] 0.57887#AllPredictions = c(PredictTest1, PredictTest2, PredictTest3)
#AllOutcomes = c(stocksTest1$PositiveDec, stocksTest2$PositiveDec, stocksTest3$PositiveDec)
#table(AllOutcomes, AllPredictions>0.5)What is the overall test-set accuracy of the cluster-then-predict approach, again using a threshold of 0.5?
We see a modest improvement over the original logistic regression model. Since predicting stock returns is a notoriously hard problem, this is a good increase in accuracy. By investing in stocks for which we are more confident that they will have positive returns (by selecting the ones with higher predicted probabilities), this cluster-then-predict model can give us an edge over the original logistic regression model.