diff --git a/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.Rmd b/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.Rmd new file mode 100644 index 00000000..e2446404 --- /dev/null +++ b/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.Rmd @@ -0,0 +1,439 @@ +--- +title: "Sessions 3-4" +author: "T. Evgeniou" +output: html_document +--- + +
+ +The purpose of this session is to become familiar with: + +1. Some visualization tools; +2. Principal Component Analysis and Factor Analysis; +3. Clustering Methods; +4. Introduction to machine learning methods; +5. A market segmentation case study. + +As always, before starting, make sure you have pulled the [session 3-4 files](https://github.com/InseadDataAnalytics/INSEADAnalytics/tree/master/CourseSessions/Sessions23) (yes, I know, it says session 2, but it is 3-4 - need to update all filenames some time, but till then we use common sense and ignore a bit the filenames) on your github repository (if you pull the course github repository you also get the session files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue). For example, assuming we are now in the "MYDIRECTORY/INSEADAnalytics" directory, we can do these: + +```{r echo=TRUE, eval=FALSE, tidy=TRUE} +getwd() +setwd("CourseSessions/Sessions23") +list.files() +rm(list=ls()) # Clean up the memory, if we want to rerun from scratch +``` +As always, you can use the `help` command in Rstudio to find out about any R function (e.g. type `help(list.files)` to learn what the R function `list.files` does). + +Let's start. + +
+
+ +### Survey Data for Market Segmentation + +We will be using the [boats case study](http://inseaddataanalytics.github.io/INSEADAnalytics/Boats-A-prerelease.pdf) as an example. At the end of this class we will be able to develop (from scratch) the readings of sessions 3-4 as well as understand the tools used and the interpretation of the results in practice - in order to make business decisions. The code used here is along the lines of the code in the session directory, e.g. in the [RunStudy.R](https://github.com/InseadDataAnalytics/INSEADAnalytics/blob/master/CourseSessions/Sessions23/RunStudy.R) file and the report [doc/Report_s23.Rmd.](https://github.com/InseadDataAnalytics/INSEADAnalytics/blob/master/CourseSessions/Sessions23/doc/Report_s23.Rmd) There may be a few differences, as there are many ways to write code to do the same thing. + +Let's load the data: + +```{r echo=FALSE, message=FALSE, prompt=FALSE, results='asis'} +source("R/library.R") +``` + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +ProjectData <- read.csv("data/Boats.csv", sep=";", dec=",") # this contains only the matrix ProjectData +ProjectData=data.matrix(ProjectData) +colnames(ProjectData)<-gsub("\\."," ",colnames(ProjectData)) +ProjectDataFactor=ProjectData[,c(2:30)] +``` +
+and do some basic visual exploration of the first 50 respondents first (always necessary to see the data first): +
+ +```{r echo=FALSE, message=FALSE, prompt=FALSE, results='asis'} +show_data = data.frame(round(ProjectData,2))[1:50,] +show_data$Variables = rownames(show_data) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` +
+ +This is the correlation matrix of the customer responses to the `r ncol(ProjectDataFactor)` attitude questions - which are the only questions that we will use for the segmentation (see the case): +
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(cbind(colnames(ProjectDataFactor), round(cor(ProjectDataFactor),2))) +m1<-gvisTable(show_data,options=list(width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE)) +print(m1,'chart') +``` +
+ +#### Questions + +1. Do you see any high correlations between the responses? Do they make sense? +2. What do these correlations imply? + +##### Answers: +
+1. I see high correlations between questions that ask for status / power / self expression +2. The correlations mean that the people interviewed will probably answer similarly to questions related to these themes + +
+
+
+ +
+ +### Key Customer Attitudes + +Clearly the survey asked many reduntant questions (can you think some reasons why?), so we may be able to actually "group" these 29 attitude questions into only a few "key factors". This not only will simplify the data, but will also greatly facilitate our understanding of the customers. + +To do so, we use methods called [Principal Component Analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) and [factor analysis](https://en.wikipedia.org/wiki/Factor_analysis) as discussed in the [session readings](http://inseaddataanalytics.github.io/INSEADAnalytics/Report_s23.html). We can use two different R commands for this (they make slightly different information easily available as output): the command `principal` (check `help(principal)` from R package [psych](http://personality-project.org/r/psych/)), and the command `PCA` from R package [FactoMineR](http://factominer.free.fr) - there are more packages and commands for this, as these methods are very widely used. + +Here is how the `principal` function is used: +
+```{r echo=TRUE, eval=TRUE, tidy=TRUE} +UnRotated_Results<-principal(ProjectDataFactor, nfactors=ncol(ProjectDataFactor), rotate="none",score=TRUE) +UnRotated_Factors<-round(UnRotated_Results$loadings,2) +UnRotated_Factors<-as.data.frame(unclass(UnRotated_Factors)) +colnames(UnRotated_Factors)<-paste("Component",1:ncol(UnRotated_Factors),sep=" ") +``` + +
+
+ +Here is how we use `PCA` one is used: +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +Variance_Explained_Table_results<-PCA(ProjectDataFactor, graph=FALSE) +Variance_Explained_Table<-Variance_Explained_Table_results$eig +Variance_Explained_Table_copy<-Variance_Explained_Table +row=1:nrow(Variance_Explained_Table) +name<-paste("Component No:",row,sep="") +Variance_Explained_Table<-cbind(name,Variance_Explained_Table) +Variance_Explained_Table<-as.data.frame(Variance_Explained_Table) +colnames(Variance_Explained_Table)<-c("Components", "Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") + +eigenvalues <- Variance_Explained_Table[,2] +``` + +
+Let's look at the **variance explained** as well as the **eigenvalues** (see session readings): +
+
+ +```{r echo=FALSE, comment=NA, warning=FALSE, error=FALSE,message=FALSE,results='asis'} +show_data = Variance_Explained_Table +m<-gvisTable(Variance_Explained_Table,options=list(width=1200, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m,'chart') +``` +
+ +```{r Fig1, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df <- cbind(as.data.frame(eigenvalues), c(1:length(eigenvalues)), rep(1, length(eigenvalues))) +colnames(df) <- c("eigenvalues", "components", "abline") +Line <- gvisLineChart(as.data.frame(df), xvar="components", yvar=c("eigenvalues","abline"), options=list(title='Scree plot', legend="right", width=900, height=600, hAxis="{title:'Number of Components', titleTextStyle:{color:'black'}}", vAxes="[{title:'Eigenvalues'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line, 'chart') +``` +
+ +#### Questions: + +1. Can you explain what this table and the plot are? What do they indicate? What can we learn from these? +2. Why does the plot have this specific shape? Could the plotted line be increasing? +3. What characteristics of these results would we prefer to see? Why? + +**Your Answers here:** +
+1. the table and the plot show the factors that explain most of the behaviour of the clients. Each factor is a combination of the questions of the survey. From the table and the plot we can learn how many factors are needed to explain a certain percentage of the behaviour (for example, to explain at least the 50% we need 5 factors). +2. the plot has this shape because the first factor is always the one that exaplain a greater part of the behaviour. +3. the smaller the number of factors, the better for the analysis of the customers. +
+
+
+ +#### Visualization and Interpretation + +Let's now see how the "top factors" look like. +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Choose one of these options: +factors_selected = sum(Variance_Explained_Table_copy[,1] >= 1) +# minimum_variance_explained = 0.5; factors_selected = 1:head(which(Variance_Explained_Table_copy[,"cumulative percentage of variance"]>= minimum_variance_explained),1) +#factors_selected = 10 +``` +
+ +To better visualise them, we will use what is called a "rotation". There are many rotations methods, we use what is called the [varimax](http://stats.stackexchange.com/questions/612/is-pca-followed-by-a-rotation-such-as-varimax-still-pca) rotation: +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Please ENTER the rotation eventually used (e.g. "none", "varimax", "quatimax", "promax", "oblimin", "simplimax", and "cluster" - see help(principal)). Defauls is "varimax" +rotation_used="varimax" +``` + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +Rotated_Results<-principal(ProjectDataFactor, nfactors=max(factors_selected), rotate=rotation_used,score=TRUE) +Rotated_Factors<-round(Rotated_Results$loadings,2) +Rotated_Factors<-as.data.frame(unclass(Rotated_Factors)) +colnames(Rotated_Factors)<-paste("Component",1:ncol(Rotated_Factors),sep=" ") +sorted_rows <- sort(Rotated_Factors[,1], decreasing = TRUE, index.return = TRUE)$ix +Rotated_Factors <- Rotated_Factors[sorted_rows,] +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis'} +show_data <- Rotated_Factors +show_data$Variables <- rownames(show_data) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` +

+ +To better visualize and interpret the factors we often "supress" loadings with small values, e.g. with absolute values smaller than 0.5. In this case our factors look as follows after suppressing the small numbers: +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +MIN_VALUE = 0.5 +Rotated_Factors_thres <- Rotated_Factors +Rotated_Factors_thres[abs(Rotated_Factors_thres) < MIN_VALUE]<-NA +colnames(Rotated_Factors_thres)<- colnames(Rotated_Factors) +rownames(Rotated_Factors_thres)<- rownames(Rotated_Factors) +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis'} +show_data <- Rotated_Factors_thres +#show_data = show_data[1:min(max_data_report,nrow(show_data)),] +show_data$Variables <- rownames(show_data) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` +

+ + +#### Questions + +1. What do the first couple of factors mean? Do they make business sense? +2. How many factors should we choose for this data/customer base? Please try a few and explain your final choice based on a) statistical arguments, b) on interpretation arguments, c) on business arguments (**you need to consider all three types of arguments**) +3. How would you interpret the factors you selected? +4. What lessons about data science do you learn when doing this analysis? Please comment. +5. (Extra/Optional) Can you make this report "dynamic" using shiny and then post it on [shinyapps.io](http://www.shinyapps.io)? (see for example exercise set 1 and interactive exercise set 2) + +**Your Answers here:** +
+1. The first factors are the ones that explain most of the customers' behaviour (for example, the first factor seems more related to the status, lifestyle, power, and the second one to feeling adventurous, socializing). They make sense because we can use them to segment and target the market. + +2. as explained before,based on eigenvalue and cumulative explained vairance, we should choose 5 factors, because they explain at least 50% of the behaviour. Besides, they don't overlap with each other, which means that with each factor we will be able to tackle one aspect of the customers preference + +3. the first factor seems more related to the status, lifestyle, power, and the second one to feeling adventurous, socializing... + +4. +
+
+
+ +
+
+ +### Market Segmentation + +Let's now use one representative question for each factor (we can also use the "factor scores" for each respondent - see [session readings](http://inseaddataanalytics.github.io/INSEADAnalytics/Report_s23.html)) to represent our survey respondents. We can choose the question with the highest absolute factor loading for each factor. For example, when we use 5 factors with the varimax rotation we can select questions Q.1.9 (I see my boat as a status symbol), Q1.18 (Boating gives me a feeling of adventure), Q1.4 (I only consider buying a boat from a reputable brand), Q1.11 (I tend to perform minor boat repairs and maintenance on my own) and Q1.2 (When buying a boat getting the lowest price is more important than the boat brand) - try it. These are columns 10, 19, 5, 12, and 3, respectively of the data matrix `Projectdata`. + +In market segmentation one may use variables to **profile** the segments which are not the same (necessarily) as those used to **segment** the market: the latter may be, for example, attitude/needs related (you define segments based on what the customers "need"), while the former may be any information that allows a company to identify the defined customer segments (e.g. demographics, location, etc). Of course deciding which variables to use for segmentation and which to use for profiling (and then **activation** of the segmentation for business purposes) is largely subjective. So in this case we will use all survey questions for profiling for now: + +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +segmentation_attributes_used = c(10,19,5,12,3) +profile_attributes_used = 2:ncol(ProjectData) +ProjectData_segment=ProjectData[,segmentation_attributes_used] +ProjectData_profile=ProjectData[,profile_attributes_used] +``` + +A key family of methods used for segmenation is what is called **clustering methods**. This is a very important problem in statistics and **machine learning**, used in all sorts of applications such as in [Amazon's pioneer work on recommender systems](http://www.cs.umd.edu/~samir/498/Amazon-Recommendations.pdf). There are many *mathematical methods* for clustering. We will use two very standard methods, **hierarchical clustering** and **k-means**. While the "math" behind all these methods can be complex, the R functions used are relatively simple to use, as we will see. + +For example, to use hierarchical clustering we simply first define some parameters used (see session readings) and then simply call the command `hclust`: + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Please ENTER the distance metric eventually used for the clustering in case of hierarchical clustering +# (e.g. "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski" - see help(dist)). +# DEFAULT is "euclidean" +distance_used="euclidean" +# Please ENTER the hierarchical clustering method to use (options are: +# "ward", "single", "complete", "average", "mcquitty", "median" or "centroid") +# DEFAULT is "ward.D" +hclust_method = "ward.D" +# Define the number of clusters: +numb_clusters_used = 3 +``` + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +Hierarchical_Cluster_distances <- dist(ProjectData_segment, method=distance_used) +Hierarchical_Cluster <- hclust(Hierarchical_Cluster_distances, method=hclust_method) + +# Assign observations (e.g. people) in their clusters +cluster_memberships_hclust <- as.vector(cutree(Hierarchical_Cluster, k=numb_clusters_used)) +cluster_ids_hclust=unique(cluster_memberships_hclust) +ProjectData_with_hclust_membership <- cbind(1:length(cluster_memberships_hclust),cluster_memberships_hclust) +colnames(ProjectData_with_hclust_membership)<-c("Observation Number","Cluster_Membership") +``` + +Finally, we can see the **dendrogram** (see class readings and online resources for more information) to have a first rough idea of what segments (clusters) we may have - and how many. +
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, fig.align='center', results='asis'} +# Display dendogram +plot(Hierarchical_Cluster, main = NULL, sub=NULL, labels = 1:nrow(ProjectData_segment), xlab="Our Observations", cex.lab=1, cex.axis=1) +# Draw dendogram with red borders around the 3 clusters +rect.hclust(Hierarchical_Cluster, k=numb_clusters_used, border="red") +``` +
+ We can also plot the "distances" traveled before we need to merge any of the lower and smaller in size clusters into larger ones - the heights of the tree branches that link the clusters as we traverse the tree from its leaves to its root. If we have n observations, this plot has n-1 numbers. +
+ + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, fig.align='center', results='asis'} +df1 <- cbind(as.data.frame(Hierarchical_Cluster$height[length(Hierarchical_Cluster$height):1]), c(1:(nrow(ProjectData)-1))) +colnames(df1) <- c("distances","index") +Line <- gvisLineChart(as.data.frame(df1), xvar="index", yvar="distances", options=list(title='Distances plot', legend="right", width=900, height=600, hAxis="{title:'Number of Components', titleTextStyle:{color:'black'}}", vAxes="[{title:'Distances'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line,'chart') +``` +
+ +To use k-means on the other hand one needs to define a priori the number of segments (which of course one can change and re-cluster). K-means also requires the choice of a few more parameters, but this is beyond our scope for now. Here is how to run K-means: +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Please ENTER the kmeans clustering method to use (options are: +# "Hartigan-Wong", "Lloyd", "Forgy", "MacQueen" +# DEFAULT is "Lloyd" +kmeans_method = "Lloyd" +# Define the number of clusters: +numb_clusters_used = 3 +kmeans_clusters <- kmeans(ProjectData_segment,centers= numb_clusters_used, iter.max=2000, algorithm=kmeans_method) +ProjectData_with_kmeans_membership <- cbind(1:length(kmeans_clusters$cluster),kmeans_clusters$cluster) +colnames(ProjectData_with_kmeans_membership)<-c("Observation Number","Cluster_Membership") + +# Assign observations (e.g. people) in their clusters +cluster_memberships_kmeans <- kmeans_clusters$cluster +cluster_ids_kmeans <- unique(cluster_memberships_kmeans) +``` + +K-means does not provide much information about segmentation. However, when we profile the segments we can start getting a better (business) understanding of what is happening. **Profiling** is a central part of segmentation: this is where we really get to mix technical and business creativity. + + +### Profiling + +There are many ways to do the profiling of the segments. For example, here we show how the *average* answers of the respondents *in each segment* compare to the *average answer of all respondents* using the ratio of the two. The idea is that if in a segment the average response to a question is very different (e.g. away from ratio of 1) than the overall average, then that question may indicate something about the segment relative to the total population. + +Here are for example the profiles of the segments using the clusters found above: + +
+ First let's see just the average answer people gave to each question for the different segments as well as the total population: +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Select whether to use the Hhierarchical clustering or the k-means clusters: + +cluster_memberships <- cluster_memberships_hclust +cluster_ids <- cluster_ids_hclust +# here is the k-means: uncomment these 2 lines +#cluster_memberships <- cluster_memberships_kmeans +#cluster_ids <- cluster_ids_kmeans + +population_average = matrix(apply(ProjectData_profile, 2, mean), ncol=1) +colnames(population_average) <- "Population" +Cluster_Profile_mean <- sapply(sort(cluster_ids), function(i) apply(ProjectData_profile[(cluster_memberships==i), ], 2, mean)) +if (ncol(ProjectData_profile) <2) + Cluster_Profile_mean=t(Cluster_Profile_mean) +colnames(Cluster_Profile_mean) <- paste("Segment", 1:length(cluster_ids), sep=" ") +cluster.profile <- cbind(population_average,Cluster_Profile_mean) +``` + + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(round(cluster.profile,2)) +#show_data = show_data[1:min(max_data_report,nrow(show_data)),] +row<-rownames(show_data) +dfnew<-cbind(row,show_data) +change<-colnames(dfnew) +change[1]<-"Variables" +colnames (dfnew)<-change +m1<-gvisTable(dfnew,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') + +``` +
+ +Let's now see the relative ratios, which we can also save in a .csv and explore if (absolutely) necessary - e.g. for collaboration with people using other tools. + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +ratio_limit = 0.1 +``` +Let's see only ratios that are larger or smaller than 1 by, say, at least `r ratio_limit`. +
+ +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +population_average_matrix <- population_average[,"Population",drop=F] %*% matrix(rep(1,ncol(Cluster_Profile_mean)),nrow=1) +cluster_profile_ratios <- (ifelse(population_average_matrix==0, 0,Cluster_Profile_mean/population_average_matrix)) +colnames(cluster_profile_ratios) <- paste("Segment", 1:ncol(cluster_profile_ratios), sep=" ") +rownames(cluster_profile_ratios) <- colnames(ProjectData)[profile_attributes_used] +## printing the result in a clean-slate table +``` + +```{r echo=TRUE, eval=TRUE, tidy=TRUE} +# Save the segment profiles in a file: enter the name of the file! +profile_file = "my_segmentation_profiles.csv" +write.csv(cluster_profile_ratios,file=profile_file) +# We can also save the cluster membership of our respondents: +data_with_segment_membership = cbind(cluster_memberships,ProjectData) +colnames(data_with_segment_membership)[1] = "Segment" +cluster_file = "my_segments.csv" +write.csv(data_with_segment_membership,file=cluster_file) +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +#library(shiny) # need this library for heatmaps to work! +# Please enter the minimum distance from "1" the profiling values should have in order to be colored +# (e.g. using heatmin = 0 will color everything - try it) +#heatmin = 0.1 +#source("R/heatmapOutput.R") +#cat(renderHeatmapX(cluster_profile_ratios, border=1, center = 1, minvalue = heatmin)) +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis'} +cluster_profile_ratios[abs(cluster_profile_ratios-1) < ratio_limit] <- NA +show_data = data.frame(round(cluster_profile_ratios,2)) +show_data$Variables <- rownames(show_data) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` + +
+
+**The further a ratio is from 1, the more important that attribute is for a segment relative to the total population.** + +
+ +#### Questions + +1. How many segments are there in our market? Why you chose that number of segments? Again, try a few and explain your final choice based on a) statistical arguments, b) on interpretation arguments, c) on business arguments (**you need to consider all three types of arguments**) +2. Can you describe the segments you found based on the profiles? +3. What if you change the number of factors and in general you *iterate the whole analysis*? **Iterations** are key in data science. +4. Can you now answer the [Boats case questions](http://inseaddataanalytics.github.io/INSEADAnalytics/Boats-A-prerelease.pdf)? What business decisions do you recommend to this company based on your analysis? + +
+ +**Your Answers here:** +
+
+
+
+ +
+ +**You have now completed your first market segmentation project.** Do you have data from another survey you can use with this report now? + +**Extra question**: explore and report a new segmentation analysis... + +... and as always **Have Fun** \ No newline at end of file diff --git a/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.html b/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.html new file mode 100644 index 00000000..acb3c35e --- /dev/null +++ b/CourseSessions/Sessions23/CopyOfSession2inclass_Ana.html @@ -0,0 +1,19472 @@ + + + + + + + + + + + + + +Sessions 3-4 + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +


+

The purpose of this session is to become familiar with:

+
    +
  1. Some visualization tools;
  2. +
  3. Principal Component Analysis and Factor Analysis;
  4. +
  5. Clustering Methods;
  6. +
  7. Introduction to machine learning methods;
  8. +
  9. A market segmentation case study.
  10. +
+

As always, before starting, make sure you have pulled the session 3-4 files (yes, I know, it says session 2, but it is 3-4 - need to update all filenames some time, but till then we use common sense and ignore a bit the filenames) on your github repository (if you pull the course github repository you also get the session files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue). For example, assuming we are now in the “MYDIRECTORY/INSEADAnalytics” directory, we can do these:

+
getwd()
+setwd("CourseSessions/Sessions23")
+list.files()
+rm(list = ls())  # Clean up the memory, if we want to rerun from scratch
+

As always, you can use the help command in Rstudio to find out about any R function (e.g. type help(list.files) to learn what the R function list.files does).

+

Let’s start.

+
+
+
+

Survey Data for Market Segmentation

+

We will be using the boats case study as an example. At the end of this class we will be able to develop (from scratch) the readings of sessions 3-4 as well as understand the tools used and the interpretation of the results in practice - in order to make business decisions. The code used here is along the lines of the code in the session directory, e.g. in the RunStudy.R file and the report doc/Report_s23.Rmd. There may be a few differences, as there are many ways to write code to do the same thing.

+

Let’s load the data:

+
ProjectData <- read.csv("data/Boats.csv", sep = ";", dec = ",")  # this contains only the matrix ProjectData
+ProjectData = data.matrix(ProjectData)
+colnames(ProjectData) <- gsub("\\.", " ", colnames(ProjectData))
+ProjectDataFactor = ProjectData[, c(2:30)]
+


and do some basic visual exploration of the first 50 respondents first (always necessary to see the data first):

+ + + + + + + +
+ +
+


+

This is the correlation matrix of the customer responses to the 29 attitude questions - which are the only questions that we will use for the segmentation (see the case):

+ + + + + + + +
+ +
+


+
+

Questions

+
    +
  1. Do you see any high correlations between the responses? Do they make sense?
  2. +
  3. What do these correlations imply?
  4. +
+
+
Answers:
+


1. I see high correlations between questions that ask for status / power / self expression 2. The correlations mean that the people interviewed will probably answer similarly to questions related to these themes

+




+
+
+
+
+
+

Key Customer Attitudes

+

Clearly the survey asked many reduntant questions (can you think some reasons why?), so we may be able to actually “group” these 29 attitude questions into only a few “key factors”. This not only will simplify the data, but will also greatly facilitate our understanding of the customers.

+

To do so, we use methods called Principal Component Analysis and factor analysis as discussed in the session readings. We can use two different R commands for this (they make slightly different information easily available as output): the command principal (check help(principal) from R package psych), and the command PCA from R package FactoMineR - there are more packages and commands for this, as these methods are very widely used.

+

Here is how the principal function is used:

+
UnRotated_Results <- principal(ProjectDataFactor, nfactors = ncol(ProjectDataFactor), 
+    rotate = "none", score = TRUE)
+UnRotated_Factors <- round(UnRotated_Results$loadings, 2)
+UnRotated_Factors <- as.data.frame(unclass(UnRotated_Factors))
+colnames(UnRotated_Factors) <- paste("Component", 1:ncol(UnRotated_Factors), 
+    sep = " ")
+



+

Here is how we use PCA one is used:

+
Variance_Explained_Table_results <- PCA(ProjectDataFactor, graph = FALSE)
+Variance_Explained_Table <- Variance_Explained_Table_results$eig
+Variance_Explained_Table_copy <- Variance_Explained_Table
+row = 1:nrow(Variance_Explained_Table)
+name <- paste("Component No:", row, sep = "")
+Variance_Explained_Table <- cbind(name, Variance_Explained_Table)
+Variance_Explained_Table <- as.data.frame(Variance_Explained_Table)
+colnames(Variance_Explained_Table) <- c("Components", "Eigenvalue", "Percentage_of_explained_variance", 
+    "Cumulative_percentage_of_explained_variance")
+
+eigenvalues <- Variance_Explained_Table[, 2]
+


Let’s look at the variance explained as well as the eigenvalues (see session readings):

+ + + + + + + +
+ +
+


+ + + + + + + +
+ +
+


+
+

Questions:

+
    +
  1. Can you explain what this table and the plot are? What do they indicate? What can we learn from these?
  2. +
  3. Why does the plot have this specific shape? Could the plotted line be increasing?
  4. +
  5. What characteristics of these results would we prefer to see? Why?
  6. +
+

Your Answers here:
1. the table and the plot show the factors that explain most of the behaviour of the clients. Each factor is a combination of the questions of the survey. From the table and the plot we can learn how many factors are needed to explain a certain percentage of the behaviour (for example, to explain at least the 50% we need 5 factors). 2. the plot has this shape because the first factor is always the one that exaplain a greater part of the behaviour. 3. the smaller the number of factors, the better for the analysis of the customers.


+
+
+

Visualization and Interpretation

+

Let’s now see how the “top factors” look like.

+
# Choose one of these options:
+factors_selected = sum(Variance_Explained_Table_copy[, 1] >= 1)
+# minimum_variance_explained = 0.5; factors_selected =
+# 1:head(which(Variance_Explained_Table_copy[,'cumulative percentage of
+# variance']>= minimum_variance_explained),1) factors_selected = 10
+


+

To better visualise them, we will use what is called a “rotation”. There are many rotations methods, we use what is called the varimax rotation:

+
# Please ENTER the rotation eventually used (e.g. 'none', 'varimax',
+# 'quatimax', 'promax', 'oblimin', 'simplimax', and 'cluster' - see
+# help(principal)). Defauls is 'varimax'
+rotation_used = "varimax"
+
Rotated_Results <- principal(ProjectDataFactor, nfactors = max(factors_selected), 
+    rotate = rotation_used, score = TRUE)
+Rotated_Factors <- round(Rotated_Results$loadings, 2)
+Rotated_Factors <- as.data.frame(unclass(Rotated_Factors))
+colnames(Rotated_Factors) <- paste("Component", 1:ncol(Rotated_Factors), sep = " ")
+sorted_rows <- sort(Rotated_Factors[, 1], decreasing = TRUE, index.return = TRUE)$ix
+Rotated_Factors <- Rotated_Factors[sorted_rows, ]
+ + + + + + + +
+ +
+



+

To better visualize and interpret the factors we often “supress” loadings with small values, e.g. with absolute values smaller than 0.5. In this case our factors look as follows after suppressing the small numbers:

+
MIN_VALUE = 0.5
+Rotated_Factors_thres <- Rotated_Factors
+Rotated_Factors_thres[abs(Rotated_Factors_thres) < MIN_VALUE] <- NA
+colnames(Rotated_Factors_thres) <- colnames(Rotated_Factors)
+rownames(Rotated_Factors_thres) <- rownames(Rotated_Factors)
+ + + + + + + +
+ +
+



+
+
+

Questions

+
    +
  1. What do the first couple of factors mean? Do they make business sense?
  2. +
  3. How many factors should we choose for this data/customer base? Please try a few and explain your final choice based on a) statistical arguments, b) on interpretation arguments, c) on business arguments (you need to consider all three types of arguments)
  4. +
  5. How would you interpret the factors you selected?
  6. +
  7. What lessons about data science do you learn when doing this analysis? Please comment.
  8. +
  9. (Extra/Optional) Can you make this report “dynamic” using shiny and then post it on shinyapps.io? (see for example exercise set 1 and interactive exercise set 2)
  10. +
+

Your Answers here:
1. The first factors are the ones that explain most of the customers’ behaviour (for example, the first factor seems more related to the status, lifestyle, power, and the second one to feeling adventurous, socializing). They make sense because we can use them to segment and target the market.

+
    +
  1. as explained before,based on eigenvalue and cumulative explained vairance, we should choose 5 factors, because they explain at least 50% of the behaviour. Besides, they don’t overlap with each other, which means that with each factor we will be able to tackle one aspect of the customers preference

  2. +
  3. the first factor seems more related to the status, lifestyle, power, and the second one to feeling adventurous, socializing…

  4. +



  5. +
+
+
+
+
+
+

Market Segmentation

+

Let’s now use one representative question for each factor (we can also use the “factor scores” for each respondent - see session readings) to represent our survey respondents. We can choose the question with the highest absolute factor loading for each factor. For example, when we use 5 factors with the varimax rotation we can select questions Q.1.9 (I see my boat as a status symbol), Q1.18 (Boating gives me a feeling of adventure), Q1.4 (I only consider buying a boat from a reputable brand), Q1.11 (I tend to perform minor boat repairs and maintenance on my own) and Q1.2 (When buying a boat getting the lowest price is more important than the boat brand) - try it. These are columns 10, 19, 5, 12, and 3, respectively of the data matrix Projectdata.

+

In market segmentation one may use variables to profile the segments which are not the same (necessarily) as those used to segment the market: the latter may be, for example, attitude/needs related (you define segments based on what the customers “need”), while the former may be any information that allows a company to identify the defined customer segments (e.g. demographics, location, etc). Of course deciding which variables to use for segmentation and which to use for profiling (and then activation of the segmentation for business purposes) is largely subjective. So in this case we will use all survey questions for profiling for now:

+


+
segmentation_attributes_used = c(10, 19, 5, 12, 3)
+profile_attributes_used = 2:ncol(ProjectData)
+ProjectData_segment = ProjectData[, segmentation_attributes_used]
+ProjectData_profile = ProjectData[, profile_attributes_used]
+

A key family of methods used for segmenation is what is called clustering methods. This is a very important problem in statistics and machine learning, used in all sorts of applications such as in Amazon’s pioneer work on recommender systems. There are many mathematical methods for clustering. We will use two very standard methods, hierarchical clustering and k-means. While the “math” behind all these methods can be complex, the R functions used are relatively simple to use, as we will see.

+

For example, to use hierarchical clustering we simply first define some parameters used (see session readings) and then simply call the command hclust:

+
# Please ENTER the distance metric eventually used for the clustering in
+# case of hierarchical clustering (e.g. 'euclidean', 'maximum', 'manhattan',
+# 'canberra', 'binary' or 'minkowski' - see help(dist)).  DEFAULT is
+# 'euclidean'
+distance_used = "euclidean"
+# Please ENTER the hierarchical clustering method to use (options are:
+# 'ward', 'single', 'complete', 'average', 'mcquitty', 'median' or
+# 'centroid') DEFAULT is 'ward.D'
+hclust_method = "ward.D"
+# Define the number of clusters:
+numb_clusters_used = 3
+
Hierarchical_Cluster_distances <- dist(ProjectData_segment, method = distance_used)
+Hierarchical_Cluster <- hclust(Hierarchical_Cluster_distances, method = hclust_method)
+
+# Assign observations (e.g. people) in their clusters
+cluster_memberships_hclust <- as.vector(cutree(Hierarchical_Cluster, k = numb_clusters_used))
+cluster_ids_hclust = unique(cluster_memberships_hclust)
+ProjectData_with_hclust_membership <- cbind(1:length(cluster_memberships_hclust), 
+    cluster_memberships_hclust)
+colnames(ProjectData_with_hclust_membership) <- c("Observation Number", "Cluster_Membership")
+

Finally, we can see the dendrogram (see class readings and online resources for more information) to have a first rough idea of what segments (clusters) we may have - and how many.

+


We can also plot the “distances” traveled before we need to merge any of the lower and smaller in size clusters into larger ones - the heights of the tree branches that link the clusters as we traverse the tree from its leaves to its root. If we have n observations, this plot has n-1 numbers.

+ + + + + + + +
+ +
+


+

To use k-means on the other hand one needs to define a priori the number of segments (which of course one can change and re-cluster). K-means also requires the choice of a few more parameters, but this is beyond our scope for now. Here is how to run K-means:

+
# Please ENTER the kmeans clustering method to use (options are:
+# 'Hartigan-Wong', 'Lloyd', 'Forgy', 'MacQueen' DEFAULT is 'Lloyd'
+kmeans_method = "Lloyd"
+# Define the number of clusters:
+numb_clusters_used = 3
+kmeans_clusters <- kmeans(ProjectData_segment, centers = numb_clusters_used, 
+    iter.max = 2000, algorithm = kmeans_method)
+ProjectData_with_kmeans_membership <- cbind(1:length(kmeans_clusters$cluster), 
+    kmeans_clusters$cluster)
+colnames(ProjectData_with_kmeans_membership) <- c("Observation Number", "Cluster_Membership")
+
+# Assign observations (e.g. people) in their clusters
+cluster_memberships_kmeans <- kmeans_clusters$cluster
+cluster_ids_kmeans <- unique(cluster_memberships_kmeans)
+

K-means does not provide much information about segmentation. However, when we profile the segments we can start getting a better (business) understanding of what is happening. Profiling is a central part of segmentation: this is where we really get to mix technical and business creativity.

+
+
+

Profiling

+

There are many ways to do the profiling of the segments. For example, here we show how the average answers of the respondents in each segment compare to the average answer of all respondents using the ratio of the two. The idea is that if in a segment the average response to a question is very different (e.g. away from ratio of 1) than the overall average, then that question may indicate something about the segment relative to the total population.

+

Here are for example the profiles of the segments using the clusters found above:

+


First let’s see just the average answer people gave to each question for the different segments as well as the total population:

+
# Select whether to use the Hhierarchical clustering or the k-means
+# clusters:
+
+cluster_memberships <- cluster_memberships_hclust
+cluster_ids <- cluster_ids_hclust
+# here is the k-means: uncomment these 2 lines cluster_memberships <-
+# cluster_memberships_kmeans cluster_ids <- cluster_ids_kmeans
+
+population_average = matrix(apply(ProjectData_profile, 2, mean), ncol = 1)
+colnames(population_average) <- "Population"
+Cluster_Profile_mean <- sapply(sort(cluster_ids), function(i) apply(ProjectData_profile[(cluster_memberships == 
+    i), ], 2, mean))
+if (ncol(ProjectData_profile) < 2) Cluster_Profile_mean = t(Cluster_Profile_mean)
+colnames(Cluster_Profile_mean) <- paste("Segment", 1:length(cluster_ids), sep = " ")
+cluster.profile <- cbind(population_average, Cluster_Profile_mean)
+ + + + + + + +
+ +
+


+

Let’s now see the relative ratios, which we can also save in a .csv and explore if (absolutely) necessary - e.g. for collaboration with people using other tools.

+
ratio_limit = 0.1
+

Let’s see only ratios that are larger or smaller than 1 by, say, at least 0.1.

+
population_average_matrix <- population_average[, "Population", drop = F] %*% 
+    matrix(rep(1, ncol(Cluster_Profile_mean)), nrow = 1)
+cluster_profile_ratios <- (ifelse(population_average_matrix == 0, 0, Cluster_Profile_mean/population_average_matrix))
+colnames(cluster_profile_ratios) <- paste("Segment", 1:ncol(cluster_profile_ratios), 
+    sep = " ")
+rownames(cluster_profile_ratios) <- colnames(ProjectData)[profile_attributes_used]
+## printing the result in a clean-slate table
+
# Save the segment profiles in a file: enter the name of the file!
+profile_file = "my_segmentation_profiles.csv"
+write.csv(cluster_profile_ratios, file = profile_file)
+# We can also save the cluster membership of our respondents:
+data_with_segment_membership = cbind(cluster_memberships, ProjectData)
+colnames(data_with_segment_membership)[1] = "Segment"
+cluster_file = "my_segments.csv"
+write.csv(data_with_segment_membership, file = cluster_file)
+ + + + + + + +
+ +
+



The further a ratio is from 1, the more important that attribute is for a segment relative to the total population.

+


+
+

Questions

+
    +
  1. How many segments are there in our market? Why you chose that number of segments? Again, try a few and explain your final choice based on a) statistical arguments, b) on interpretation arguments, c) on business arguments (you need to consider all three types of arguments)
  2. +
  3. Can you describe the segments you found based on the profiles?
  4. +
  5. What if you change the number of factors and in general you iterate the whole analysis? Iterations are key in data science.
  6. +
  7. Can you now answer the Boats case questions? What business decisions do you recommend to this company based on your analysis?
  8. +
+


+

Your Answers here:



+
+

You have now completed your first market segmentation project. Do you have data from another survey you can use with this report now?

+

Extra question: explore and report a new segmentation analysis…

+

… and as always Have Fun

+
+
+ + +
+ + + + + + + + diff --git a/Exercises/Exerciseset1/DataSet1.Rdata b/Exercises/Exerciseset1/DataSet1.Rdata new file mode 100644 index 00000000..48f5ae6b Binary files /dev/null and b/Exercises/Exerciseset1/DataSet1.Rdata differ diff --git a/Exercises/Exerciseset1/ExerciseSet1.Rmd b/Exercises/Exerciseset1/ExerciseSet1.Rmd index 93d3f040..4aa35892 100644 --- a/Exercises/Exerciseset1/ExerciseSet1.Rmd +++ b/Exercises/Exerciseset1/ExerciseSet1.Rmd @@ -1,7 +1,11 @@ --- -title: "Exercise Set 1" +title: "Exercise Set 1 - Solutions Ana de Andres" author: "T. Evgeniou" +<<<<<<< HEAD +date: "7 Jan 2016" +======= +>>>>>>> refs/remotes/InseadDataAnalytics/master output: html_document --- @@ -85,6 +89,14 @@ pnl_plot(SPY) 3. *(Extra points)* What if we want to also see the returns of another company, say Yahoo!, in the same period? Can you get that data and report the statistics for Yahoo!'s stock, too? **Your Answers here:** + +1. in dataSet1.R, the variable mytickers = c("SPY", "AAPL") - line 9 - includes AAPL, which is the trading name of Apple. THe lines from 13 to 46 download the data for the values on mytickers, which include the whole market and Apple specifically. + +2. For this question I created a variable called AAPL = StockReturns[,"AAPL"]. The cumulative return is 384.1, the average daily return is 0.139 and the std deviation is 2.188 (using the same functions as with SPY: round(100*sum(AAPL),1), round(100*mean(AAPL),3) and round(100*sd(AAPL),3)) +We can get the same results if we get the data from the second column of Stockreturns, which stores the values of Apple (example for cumulative round(100*sum(StockReturns[,2]),1)) + +3. We could get all the different companies traded from Yahoo. We would need to add the trading names to the variable mytickers. We would then store the data in additional columns in both StockPrices and StockReturns. +


@@ -126,6 +138,9 @@ We can also transpose the matrix of returns to create a new "horizontal" matrix. 2. Based on the help for the R function *apply* (`help(apply)`), can you create again the portfolio of S&P and Apple and plot the returns in a new figure below? **Your Answers here:** +1. We use the command dim(transposedData) to get all the dimensions (which returs 2 2771). If we want to get only the numer of rows we use nrow(transposedData) and for the number of columns we use ncol(transposedData), which gives us again 2 and 2771. + +2. I assume we are using the new variable transposedData for this question. With this variable, we would tell R to apply the sum over the columns, which is given by a 2 in the Margin argument. The code would then be portfolio2 = apply(transposedData, 2, mean). We can then plot the new variable portfolio2: pnl_plot(portfolio2).


@@ -144,6 +159,7 @@ This is an important step and will get you to think about the overall process on 2. *(Extra Exercise)* Can you get the returns of a few companies and plot the returns of an equal weighted portfolio with those companies during some period you select? **Your Answers here:** +1. We should change line 11 in DataSet1.R, which includes the variable startDate = "2005-01-01". We would assign then startDate = "2001-01-01"


@@ -178,6 +194,7 @@ myData + StockReturns[1:40,]) ``` **Your Answers here:** +1. It reutrns a value of 20.


diff --git a/Exercises/Exerciseset1/ExerciseSet1.html b/Exercises/Exerciseset1/ExerciseSet1.html new file mode 100644 index 00000000..fd82844a --- /dev/null +++ b/Exercises/Exerciseset1/ExerciseSet1.html @@ -0,0 +1,216 @@ + + + + + + + + + + + + + +Exercise Set 1 - Solutions Ana de Andres + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +


+

The purpose of this exercise is to become familiar with:

+
    +
  1. Basic statistics functions in R;
  2. +
  3. Simple matrix operations;
  4. +
  5. Simple data manipulations;
  6. +
  7. The idea of functions as well as some useful customized functions provided.
  8. +
+

While doing this exercise we will also see how to generate replicable and customizable reports. For this purpose the exercise uses the R Markdown capabilities (see Markdown Cheat Sheet or a basic introduction to R Markdown). These capabilities allow us to create dynamic reports. For example today’s date is 2016-01-07 (you need to see the .Rmd to understand that this is not a static typed-in date but it changes every time you compile the .Rmd - if the date changed of course).

+

Before starting, make sure you have pulled the exercise files on your github repository (if you pull the course github repository you also get the exercise set files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue):

+
getwd()
+setwd("Exercises/Exerciseset1/")
+list.files()
+

Note: you can always use the help command in Rstudio to find out about any R function (e.g. type help(list.files) to learn what the R function list.files does).

+

Let’s now see the exercise.

+

IMPORTANT: You should answer all questions by simply adding your code/answers in this document through editing the file ExerciseSet1.Rmd and then clicking on the “Knit HTML” button in RStudio. Once done, please post your .Rmd and html files in your github repository.

+
+

Exercise Data

+

We download daily prices (open, high, low, close, and adjusted close) and volume data of publicly traded companies and markets from the web (e.g. Yahoo! or Google, etc). This is done by sourcing the file data.R as well as some helper functions in herpersSet1.R which also installs a number of R libraries (hence the first time you run this code you will see a lot of red color text indicating the download and installation process):

+
source("helpersSet1.R")
+source("dataSet1.R")
+

We have 2771 days of data, starting from 2005-01-04 until 2016-01-06.

+
+
+

Part I: Statistics of S&P Daily Returns

+

Here are some basic statistics about the S&P returns:

+
    +
  1. The cumulative returns of the S&P index during this period is 94.4%.
  2. +
  3. The average daily returns of the S&P index during this period is 0.034%;
  4. +
  5. The standard deviation of the daily returns of the S&P index during this period is 1.257%;
  6. +
+

Here are returns of the S&P in this period (note the use of the helper function pnl_plot - defined in file helpersSet1.R):

+

+
+

Questions

+
    +
  1. Notice that the code also downloads the returns of Apple during the same period. Can you explain where this is done in the code (including the .R files used)?
  2. +
  3. What are the cumulative, average daily returns, and the standard deviation of the daily returns of Apple in the same period?
  4. +
  5. (Extra points) What if we want to also see the returns of another company, say Yahoo!, in the same period? Can you get that data and report the statistics for Yahoo!’s stock, too?
  6. +
+

Your Answers here:

+
    +
  1. in dataSet1.R, the variable mytickers = c(“SPY”, “AAPL”) - line 9 - includes AAPL, which is the trading name of Apple. THe lines from 13 to 46 download the data for the values on mytickers, which include the whole market and Apple specifically.

  2. +
  3. For this question I created a variable called AAPL = StockReturns[,“AAPL”]. The cumulative return is 384.1, the average daily return is 0.139 and the std deviation is 2.188 (using the same functions as with SPY: round(100sum(AAPL),1), round(100mean(AAPL),3) and round(100sd(AAPL),3)) We can get the same results if we get the data from the second column of Stockreturns, which stores the values of Apple (example for cumulative round(100sum(StockReturns[,2]),1))

  4. +
  5. We could get all the different companies traded from Yahoo. We would need to add the trading names to the variable mytickers. We would then store the data in additional columns in both StockPrices and StockReturns.

  6. +
+



+
+
+
+

Part II: Simple Matrix Manipulations

+

For this part of the exercise we will do some basic manipulations of the data. First note that the data are in a so-called matrix format. If you run these commands in RStudio (use help to find out what they do) you will see how matrices work:

+
class(StockReturns)
+dim(StockReturns)
+nrow(StockReturns)
+ncol(StockReturns)
+StockReturns[1:4,]
+head(StockReturns,5)
+tail(StockReturns,5) 
+

We will now use an R function for matrices that is extremely useful for analyzing data. It is called apply. Check it out using help in R.

+

For example, we can now quickly estimate the average returns of S&P and Apple (of course this can be done manually, too, but what if we had 500 stocks - e.g. a matrix with 500 columns?) and plot the returns of that 50-50 on S&P and Apple portfolio:

+

+

We can also transpose the matrix of returns to create a new “horizontal” matrix. Let’s call this matrix (variable name) transposedData. We can do so using this command: transposedData = t(StockReturns).

+
+

Questions

+
    +
  1. What R commands can you use to get the number of rows and number of columns of the new matrix called transposedData?
  2. +
  3. Based on the help for the R function apply (help(apply)), can you create again the portfolio of S&P and Apple and plot the returns in a new figure below?
  4. +
+

Your Answers here: 1. We use the command dim(transposedData) to get all the dimensions (which returs 2 2771). If we want to get only the numer of rows we use nrow(transposedData) and for the number of columns we use ncol(transposedData), which gives us again 2 and 2771.

+
    +
  1. I assume we are using the new variable transposedData for this question. With this variable, we would tell R to apply the sum over the columns, which is given by a 2 in the Margin argument. The code would then be portfolio2 = apply(transposedData, 2, mean). We can then plot the new variable portfolio2: pnl_plot(portfolio2).

  2. +
+
+
+
+

Part III: Reproducibility and Customization

+

This is an important step and will get you to think about the overall process once again.

+
+

Questions

+
    +
  1. We want to re-do all this analysis with data since 2001-01-01: what change do we need to make in the code (hint: all you need to change is one line - exactly 1 number! - in data.R file), and how can you get the new exercise set with the data since 2001-01-01?
  2. +
  3. (Extra Exercise) Can you get the returns of a few companies and plot the returns of an equal weighted portfolio with those companies during some period you select?
  4. +
+

Your Answers here: 1. We should change line 11 in DataSet1.R, which includes the variable startDate = “2005-01-01”. We would assign then startDate = “2001-01-01”

+
+
+
+

Part IV: Read/Write .CSV files

+

Finally, one can read and write data in .CSV files. For example, we can save the first 20 days of data for S&P and Apple in a file using the command:

+
write.csv(StockReturns[1:20,c("SPY","AAPL")], file = "twentydays.csv", row.names = TRUE, col.names = TRUE) 
+

Do not get surpsised if you see the csv file in your directories suddenly! You can then read the data from the csv file using the read.csv command. For example, this will load the data from the csv file and save it in a new variable that now is called “myData”:

+
myData <- read.csv(file = "twentydays.csv", header = TRUE, sep=";")
+

Try it!

+
+

Questions

+
    +
  1. Once you write and read the data as described above, what happens when you run this command in the console of the RStudio: sum(myData != StockReturns[1:20,])
  2. +
  3. (Extra exercise) What do you think will happen if you now run this command, and why:
  4. +
+
myData + StockReturns[1:40,])
+

Your Answers here: 1. It reutrns a value of 20.

+
+
+
+

Extra Question

+

Can you now load another dataset from some CSV file and report some basic statistics about that data?

+


+
+
+

Creating Interactive Documents

+

Finally, just for fun, one can add some interactivity in the report using Shiny.All one needs to do is set the eval flag of the code chunk below (see the .Rmd file) to “TRUE”, add the line “runtime: shiny” at the very begining of the .Rmd file, make the markdown output to be “html_document”, and then press “Run Document”.

+
sliderInput("startdate", "Starting Date:", min = 1, max = length(portfolio), 
+            value = 1)
+sliderInput("enddate", "End Date:", min = 1, max = length(portfolio), 
+            value = length(portfolio))
+
+renderPlot({
+  pnl_plot(portfolio[input$startdate:input$enddate])
+})
+

Have fun.

+
+ + +
+ + + + + + + + diff --git a/Exercises/Exerciseset1/ExerciseSet1_Ana de Andres.Rmd b/Exercises/Exerciseset1/ExerciseSet1_Ana de Andres.Rmd new file mode 100644 index 00000000..c42c0be2 --- /dev/null +++ b/Exercises/Exerciseset1/ExerciseSet1_Ana de Andres.Rmd @@ -0,0 +1,184 @@ + +--- +title: "Exercise Set 1 - Solutions Ana de Andres " +author: "T. Evgeniou" +date: "7 Jan 2016" +output: html_document +--- + + +
+ +The purpose of this exercise is to become familiar with: + +1. Basic statistics functions in R; +2. Simple matrix operations; +3. Simple data manipulations; +4. The idea of functions as well as some useful customized functions provided. + +While doing this exercise we will also see how to generate replicable and customizable reports. For this purpose the exercise uses the R Markdown capabilities (see [Markdown Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/rmarkdown-cheatsheet.pdf) or a [basic introduction to R Markdown](http://rmarkdown.rstudio.com/authoring_basics.html)). These capabilities allow us to create dynamic reports. For example today's date is `r Sys.Date()` (you need to see the .Rmd to understand that this is *not* a static typed-in date but it changes every time you compile the .Rmd - if the date changed of course). + +Before starting, make sure you have pulled the [exercise files](https://github.com/InseadDataAnalytics/INSEADAnalytics/tree/master/Exercises/Exerciseset1) on your github repository (if you pull the course github repository you also get the exercise set files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue): + +```{r echo=TRUE, eval=FALSE, tidy=TRUE} +getwd() +setwd("Exercises/Exerciseset1/") +list.files() +``` + +**Note:** you can always use the `help` command in Rstudio to find out about any R function (e.g. type `help(list.files)` to learn what the R function `list.files` does). + +Let's now see the exercise. + +**IMPORTANT:** You should answer all questions by simply adding your code/answers in this document through editing the file ExerciseSet1.Rmd and then clicking on the "Knit HTML" button in RStudio. Once done, please post your .Rmd and html files in your github repository. + +### Exercise Data + +We download daily prices (open, high, low, close, and adjusted close) and volume data of publicly traded companies and markets from the web (e.g. Yahoo! or Google, etc). This is done by sourcing the file data.R as well as some helper functions in herpersSet1.R which also installs a number of R libraries (hence the first time you run this code you will see a lot of red color text indicating the *download* and *installation* process): + +```{r eval = TRUE, echo=TRUE, error = FALSE, warning=FALSE,message=FALSE,results='asis'} +source("helpersSet1.R") +source("dataSet1.R") +``` + +We have `r nrow(StockReturns)` days of data, starting from `r rownames(StockReturns)[1]` until `r tail(rownames(StockReturns),1)`. + +### Part I: Statistics of S&P Daily Returns + +Here are some basic statistics about the S&P returns: + +1. The cumulative returns of the S&P index during this period is `r round(100*sum(StockReturns[,1]),1)`%. +2. The average daily returns of the S&P index during this period is `r round(100*mean(StockReturns[,1]),3)`%; +2. The standard deviation of the daily returns of the S&P index during this period is `r round(100*sd(StockReturns[,1]),3)`%; + +Here are returns of the S&P in this period (note the use of the helper function pnl_plot - defined in file helpersSet1.R): + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis',fig.align='center', fig.height=4,fig.width= 6, fig=TRUE} +SPY = StockReturns[,"SPY"] +pnl_plot(SPY) +``` + +#### Questions + +1. Notice that the code also downloads the returns of Apple during the same period. Can you explain where this is done in the code (including the .R files used)? +2. What are the cumulative, average daily returns, and the standard deviation of the daily returns of Apple in the same period? +3. *(Extra points)* What if we want to also see the returns of another company, say Yahoo!, in the same period? Can you get that data and report the statistics for Yahoo!'s stock, too? + +**Your Answers here:** + +1. in dataSet1.R, the variable mytickers = c("SPY", "AAPL") - line 9 - includes AAPL, which is the trading name of Apple. THe lines from 13 to 46 download the data for the values on mytickers, which include the whole market and Apple specifically. + +2. For this question I created a variable called AAPL = StockReturns[,"AAPL"]. The cumulative return is 384.1, the average daily return is 0.139 and the std deviation is 2.188 (using the same functions as with SPY: round(100*sum(AAPL),1), round(100*mean(AAPL),3) and round(100*sd(AAPL),3)) +We can get the same results if we get the data from the second column of Stockreturns, which stores the values of Apple (example for cumulative round(100*sum(StockReturns[,2]),1)) + +3. We could get all the different companies traded from Yahoo. We would need to add the trading names to the variable mytickers. We would then store the data in additional columns in both StockPrices and StockReturns. + +
+
+ +### Part II: Simple Matrix Manipulations + +For this part of the exercise we will do some basic manipulations of the data. First note that the data are in a so-called matrix format. If you run these commands in RStudio (use help to find out what they do) you will see how matrices work: + +```{r eval = FALSE, echo=TRUE} +class(StockReturns) +dim(StockReturns) +nrow(StockReturns) +ncol(StockReturns) +StockReturns[1:4,] +head(StockReturns,5) +tail(StockReturns,5) +``` + +We will now use an R function for matrices that is extremely useful for analyzing data. It is called *apply*. Check it out using help in R. + +For example, we can now quickly estimate the average returns of S&P and Apple (of course this can be done manually, too, but what if we had 500 stocks - e.g. a matrix with 500 columns?) and plot the returns of that 50-50 on S&P and Apple portfolio: + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis',fig.align='center', fig=TRUE} +portfolio = apply(StockReturns,1,mean) +names(portfolio) <- rownames(StockReturns) +pnl_plot(portfolio) +``` + + +We can also transpose the matrix of returns to create a new "horizontal" matrix. Let's call this matrix (variable name) transposedData. We can do so using this command: `transposedData = t(StockReturns)`. + +#### Questions + +1. What R commands can you use to get the number of rows and number of columns of the new matrix called transposedData? +2. Based on the help for the R function *apply* (`help(apply)`), can you create again the portfolio of S&P and Apple and plot the returns in a new figure below? + +**Your Answers here:** +1. We use the command dim(transposedData) to get all the dimensions (which returs 2 2771). If we want to get only the numer of rows we use nrow(transposedData) and for the number of columns we use ncol(transposedData), which gives us again 2 and 2771. + +2. I assume we are using the new variable transposedData for this question. With this variable, we would tell R to apply the sum over the columns, which is given by a 2 in the Margin argument. The code would then be portfolio2 = apply(transposedData, 2, mean). We can then plot the new variable portfolio2: pnl_plot(portfolio2). +
+
+ +### Part III: Reproducibility and Customization + +This is an important step and will get you to think about the overall process once again. + +#### Questions + +1. We want to re-do all this analysis with data since 2001-01-01: what change do we need to make in the code (hint: all you need to change is one line - exactly 1 number! - in data.R file), and how can you get the new exercise set with the data since 2001-01-01? +2. *(Extra Exercise)* Can you get the returns of a few companies and plot the returns of an equal weighted portfolio with those companies during some period you select? + +**Your Answers here:** +1. We should change line 11 in DataSet1.R, which includes the variable startDate = "2005-01-01". We would assign then startDate = "2001-01-01" +
+
+ +### Part IV: Read/Write .CSV files + +Finally, one can read and write data in .CSV files. For example, we can save the first 20 days of data for S&P and Apple in a file using the command: + +```{r eval = TRUE, echo=TRUE, comment=NA, warning=FALSE, message=FALSE,results='asis'} +write.csv(StockReturns[1:20,c("SPY","AAPL")], file = "twentydays.csv", row.names = TRUE, col.names = TRUE) +``` + +Do not get surpsised if you see the csv file in your directories suddenly! You can then read the data from the csv file using the read.csv command. For example, this will load the data from the csv file and save it in a new variable that now is called "myData": + +```{r eval = TRUE, echo=TRUE, comment=NA, warning=FALSE, message=FALSE,results='asis'} +myData <- read.csv(file = "twentydays.csv", header = TRUE, sep=";") +``` + +Try it! + +#### Questions + +1. Once you write and read the data as described above, what happens when you run this command in the console of the RStudio: `sum(myData != StockReturns[1:20,])` +2. *(Extra exercise)* What do you think will happen if you now run this command, and why: + +```{r eval = FALSE, echo=TRUE} +myData + StockReturns[1:40,]) +``` + +**Your Answers here:** +1. It reutrns a value of 20. +
+
+ +### Extra Question + +Can you now load another dataset from some CSV file and report some basic statistics about that data? + +
+ +### Creating Interactive Documents + +Finally, just for fun, one can add some interactivity in the report using [Shiny](http://rmarkdown.rstudio.com/authoring_shiny.html).All one needs to do is set the eval flag of the code chunk below (see the .Rmd file) to "TRUE", add the line "runtime: shiny" at the very begining of the .Rmd file, make the markdown output to be "html_document", and then press "Run Document". + +```{r, eval=FALSE, echo = TRUE} +sliderInput("startdate", "Starting Date:", min = 1, max = length(portfolio), + value = 1) +sliderInput("enddate", "End Date:", min = 1, max = length(portfolio), + value = length(portfolio)) + +renderPlot({ + pnl_plot(portfolio[input$startdate:input$enddate]) +}) +``` + +Have fun. + diff --git a/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.Rmd b/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.Rmd new file mode 100644 index 00000000..13cdbdea --- /dev/null +++ b/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.Rmd @@ -0,0 +1,347 @@ +--- +title: "Exercise Set 2: A $300 Billion Strategy" +author: "T. Evgeniou" +output: html_document +--- + +
+ +The purpose of this exercise is to become familiar with: + +1. Some time series analysis tools; +2. Correlation matrices and principal component analysis (PCA) (see [readings of sessions 3-4](http://inseaddataanalytics.github.io/INSEADAnalytics/Report_s23.html)); +3. More data manipulation and reporting tools (including Google Charts). + +As always, while doing this exercise we will also see how to generate replicable and customizable reports. For this purpose the exercise uses the R Markdown capabilities (see [Markdown Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/rmarkdown-cheatsheet.pdf) or a [basic introduction to R Markdown](http://rmarkdown.rstudio.com/authoring_basics.html)). These capabilities allow us to create dynamic reports. For example today's date is `r Sys.Date()` (you need to see the .Rmd to understand that this is *not* a static typed-in date but it changes every time you compile the .Rmd - if the date changed of course). + +Before starting, make sure you have pulled the [exercise set 2 souce code files](https://github.com/InseadDataAnalytics/INSEADAnalytics/tree/master/Exercises/Exerciseset2) on your github repository (if you pull the course github repository you also get the exercise set files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue). For example, assuming we are now in the "Data Analytics R version/INSEADAnalytics" directory, we can do these: + +```{r echo=TRUE, eval=FALSE, tidy=TRUE} +getwd() +setwd("Exercises/Exerciseset2/") +list.files() +``` + +**Note:** as always, you can use the `help` command in Rstudio to find out about any R function (e.g. type `help(list.files)` to learn what the R function `list.files` does). + +Let's now see the exercise. + +**IMPORTANT:** You should answer all questions by simply adding your code/answers in this document through editing the file ExerciseSet2.Rmd and then clicking on the "Knit HTML" button in RStudio. Once done, please post your .Rmd and html files in your github repository. + +
+ +### The Exercise: Introduction + +For this exercise we will use the Futures' daily returns to develop what is considered to be a *"classic" hedge fund trading strategy*, a **futures trend following strategy**. There is a lot written about this, so it is worth doing some online search about "futures trend following", or "Managed Futures", or "Commodity Trading Advisors (CTA)". There is about **[$300 billion](http://www.barclayhedge.com/research/indices/cta/Money_Under_Management.html)** invested on this strategy today, and is considered to be one of the **oldest hedge fund strategies**. Some example links are: + +* [A fascinating report on 2 centuries of trend following from the CFM hedge - a $6 billion fund](https://www.trendfollowing.com/whitepaper/Two_Centuries_Trend_Following.pdf) +* [Another fascinating report on 1 century of trend following investing from AQR - a $130 billion fund](https://www.aqr.com/library/aqr-publications/a-century-of-evidence-on-trend-following-investing) +* [Wikipedia on CTAs](https://en.wikipedia.org/wiki/Commodity_trading_advisor) +* [Morningstar on CTAs](http://www.morningstar.co.uk/uk/news/69379/commodity-trading-advisors-(cta)-explained.aspx) +* [A report](http://perspectives.pictet.com/wp-content/uploads/2011/01/Trading-Strategies-Final.pdf) +* [Man AHL (a leading hedge fund on CTAs - among others) - an $80 billion fund](https://www.ahl.com) + +Of course there are also many starting points for developing such a strategy (for example [this R bloggers one](http://www.r-bloggers.com/system-from-trend-following-factors/) (also on [github](https://gist.github.com/timelyportfolio/2855303)), or the [turtle traders website](http://turtletrader.com) which has many resources. + +In this exercise we will develop our own strategy from scratch. + +*Note (given today's market conditions):* **Prices of commodities, like oil or gold, can be excellent indicators of the health of the economy and of various industries, as we will also see below**. + +### Getting the Futures Data + +There are many ways to get futures data. For example, one can use the [Quandl package,](https://www.quandl.com/browse) or the [turtle traders resources,](http://turtletrader.com/hpd/) or (for INSEAD only) get data from the [INSEAD library finance data resources](http://sites.insead.edu/library/E_resources/ER_subject.cfm#Stockmarket) website. One has to pay attention on how to create continuous time series from underlying contracts with varying deliveries (e.g. see [here](https://www.quantstart.com/articles/Continuous-Futures-Contracts-for-Backtesting-Purposes) ). Using a combination of the resources above, we will use data for a number of commodities. + + +### Data description + +Let's load the data and see what we have. + +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +source("helpersSet2.R") +library(googleVis) +load("data/FuturesTrendFollowingData.Rdata") +``` + +
+We have data from `r head(rownames(futures_data),1)` to `r tail(rownames(futures_data),1)` of daily returns for the following `r ncol(futures_data)` futures: + +
+ +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE, results='asis'} +show_data = data.frame(colnames(futures_data)) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` +
+ + + +### Basic data analysis + +Let's see how these are correlated. Let's also make it look nicer (than, say, what we did in Exercise Set 1), using [Google Charts](https://code.google.com/p/google-motion-charts-with-r/wiki/GadgetExamples) (see examples online, e.g. [examples](https://cran.r-project.org/web/packages/googleVis/vignettes/googleVis_examples.html) and the [R package used used](https://cran.r-project.org/web/packages/googleVis/googleVis.pdf) ).The correlation matrix is as follows (note that the table is "dynamic": for example you can sort it based on each column by clicking on the column's header) + +
+ + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(cbind(colnames(futures_data), round(cor(futures_data),2))) +m1<-gvisTable(show_data,options=list(width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE)) +print(m1,'chart') +``` + +
+ +We see quite high correlations among some of the futures. Does it make sense? Why? Do you see some negative correlations? Do those make sense? + +Given such high correlations, we can try to see whether there are some "principal components" (see [reading on dimensionality reduction](http://inseaddataanalytics.github.io/INSEADAnalytics/Report_s23.html)). This analysis can also indicate whether all futures (the global economy!) are driven by some common "factors" (let's call them **"risk factors"**). + +
+ +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +Variance_Explained_Table_results<-PCA(futures_data, graph=FALSE) +Variance_Explained_Table<-cbind(paste("component",1:ncol(futures_data),sep=" "),Variance_Explained_Table_results$eig) +Variance_Explained_Table<-as.data.frame(Variance_Explained_Table) +colnames(Variance_Explained_Table)<-c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(Variance_Explained_Table) +m1<-gvisTable(show_data,options=list(width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1,'chart') +``` +
+ +Here is the scree plot (see Sessions 3-4 readings): +
+ +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +eigenvalues <- Variance_Explained_Table[,2] +``` + +```{r Fig1, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df <- cbind(as.data.frame(eigenvalues), c(1:length(eigenvalues)), rep(1, length(eigenvalues))) +colnames(df) <- c("eigenvalues", "components", "abline") +Line <- gvisLineChart(as.data.frame(df), xvar="components", yvar=c("eigenvalues","abline"), options=list(title='Scree plot', legend="right", width=900, height=600, hAxis="{title:'Number of Components', titleTextStyle:{color:'black'}}", vAxes="[{title:'Eigenvalues'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line, 'chart') +``` + +
+ +Let's now see how the 20 first (**rotated**) principal components look like. Let's also use the *rotated* factors (note that these are not really the "principal component", as explained in the [reading on dimensionality reduction](http://inseaddataanalytics.github.io/INSEADAnalytics/Report_s23.html)) and not show any numbers less than 0.3 in absolute value, to avoid cluttering. Note again that you can sort the table according to any column by clicking on the header of that column. +
+ +```{r echo=TRUE, comment=NA, warning=FALSE, error=FALSE,message=FALSE,results='asis',tidy=TRUE} +corused = cor(futures_data[,apply(futures_data!=0,2,sum) > 10, drop=F]) +Rotated_Results<-principal(corused, nfactors=20, rotate="varimax",score=TRUE) +Rotated_Factors<-round(Rotated_Results$loadings,2) +Rotated_Factors<-as.data.frame(unclass(Rotated_Factors)) +colnames(Rotated_Factors)<-paste("Component",1:ncol(Rotated_Factors),sep=" ") + +sorted_rows <- sort(Rotated_Factors[,1], decreasing = TRUE, index.return = TRUE)$ix +Rotated_Factors <- Rotated_Factors[sorted_rows,] +Rotated_Factors[abs(Rotated_Factors) < 0.3]<-NA +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, error=FALSE,message=FALSE,results='asis'} +show_data <- Rotated_Factors +show_data<-cbind(rownames(show_data),show_data) +colnames(show_data)<-c("Variables",colnames(Rotated_Factors)) +m1<-gvisTable(show_data,options=list(showRowNumber=TRUE,width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable')) +print(m1,'chart') +``` +
+ +#### Questions: + +1. How many principal components ("factors") do we need to explain at least 50% of the variance in this data? +2. What are the highest weights (in absolute value) of the first principal component portfolio above on the `r ncol(futures_data)` futures? +3. Can we interpret the first 10 components? How would you call these factors? +4. Can you now generate the principal components and scree plot using only: a) the pre-crisis bull market years (e.g. only using the data between November 1, 2002, and October 1, 2007)? b) the financial crisis years (e.g. only using the data between October 1, 2007 and March 1, 2009), (Hint: you can select subsets of the data using for example the command `crisis_data = futures_data[as.Date(rownames(futures_data)) > "2007-10-01" & as.Date(rownames(futures_data)) < "2009-03-01", ]) +5. Based on your analysis in question 3, please discuss any differences you observe about the futures returns during bull and bear markets. What implications may these results have? What do the results imply about how assets are correlated during bear years compared to bull years? +6. (Extra - optional) Can you create an interactive (shiny based) tool so that we can study how the "**risk factors**" change ove time? (Hint: see [Exercise set 1](https://github.com/InseadDataAnalytics/INSEADAnalytics/blob/master/Exercises/Exerciseset1/ExerciseSet1.Rmd) and online resources on [Shiny](http://rmarkdown.rstudio.com/authoring_shiny.html) such as these [Shiny lessons](http://shiny.rstudio.com/tutorial/lesson1/). Note however that you may need to pay attention to various details e.g. about how to include Google Charts in Shiny tools - so keep this extra exercise for later!). + +
+ +**Your Answers here:** + +1. We need 6 factors (cumulative percentage = 52,29%) + +2. In absolute values, the highest weights are 0,93 (5yr T-Notes US and 10yr T notes), 0,86 (treasury bonds), 0,85 (2 yr T-note US), 0,79, 0,74, ... + +3. If we see the first 10 factors, we can see that, since the overlapping of values between the factors are low, we can explain the variables with these 10 factors. + +4. +For the pre crisis data: +```{r echo=FALSE, comment=NA, warning=FALSE, error=FALSE,message=FALSE,results='asis'} +pre_crisis_data = futures_data[as.Date(rownames(futures_data)) > "2002-11-01" & as.Date(rownames(futures_data)) < "2007-10-01", ] +Variance_Explained_Table_results<-PCA(pre_crisis_data, graph=FALSE) +Variance_Explained_Table<-cbind(paste("component",1:ncol(pre_crisis_data),sep=" "),Variance_Explained_Table_results$eig) +Variance_Explained_Table<-as.data.frame(Variance_Explained_Table) +colnames(Variance_Explained_Table)<-c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") +show_data = data.frame(Variance_Explained_Table) +m1<-gvisTable(show_data,options=list(width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1,'chart') +``` +For the crisis data: + +```{r echo=FALSE, comment=NA, warning=FALSE, error=FALSE,message=FALSE,results='asis'} +crisis_data = futures_data[as.Date(rownames(futures_data)) > "2007-10-01" & as.Date(rownames(futures_data)) < "2009-03-01", ] + +Variance_Explained_Table_results<-PCA(crisis_data, graph=FALSE) +Variance_Explained_Table<-cbind(paste("component",1:ncol(crisis_data),sep=" "),Variance_Explained_Table_results$eig) +Variance_Explained_Table<-as.data.frame(Variance_Explained_Table) +colnames(Variance_Explained_Table)<-c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") +show_data = data.frame(Variance_Explained_Table) +m1<-gvisTable(show_data,options=list(width=1920, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1,'chart') +``` + +5. +During the bull (pre crisis) years we need more components to explain at least 50% of the data (again 6 factors, while during the crisis years we can explain it with 4 factors) + + + +
+ +### A Simple Futures Trend Following Strategy + +We can now develop a simple futures trend following trading strategy, as outlined in the papers in the Exercise Introduction above. There are about $300 billion invested in such strategies! Of course we cannot develop here a sophisticated product, but with some more work... + +We will do the following: + +1. Calculate a number of moving averages of different "window lengths" for each of the `r ncol(futures_data)` futures - there are [many](http://www.r-bloggers.com/stock-analysis-using-r/) so called [technical indicators](http://www.investopedia.com/active-trading/technical-indicators/) one can use. We will use the "moving average" function `ma` for this (try for example to see what this returns `ma(1:10,2)` ). +2. Add the signs (can also use the actual moving average values of course - try it!) of these moving averages (as if they "vote"), and then scale this sum across all futures so that the sum of their (of the sum across all futures!) absolute value across all futures is 1 (hence we invest $1 every day - you see why?). +3. Then invest every day in each of the `r ncol(futures_data)` an amount that is defined by the weights calculated in step 2, using however the weights calculated using data until 2 days ago (why 2 days and not 1 day?) - see the use of the helper function `shift` for this. +4. Finally see the performance of this strategy. + +Here is the code. +
+ +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +signal_used = 0*futures_data # just initialize the trading signal to be 0 +# Take many moving Average (MA) Signals and let them "vote" with their sign (+-1, e.g. long or short vote, for each signal) +MAfreq<-seq(20,250,by=30) +for (iter in 1:length(MAfreq)) + signal_used = signal_used + sign(apply(futures_data,2, function(r) ma(r,MAfreq[iter]))) +# Now make sure we invest $1 every day (so the sum of the absolute values of the weights is 1 every day) +signal_used = t(apply(signal_used,1,function(r) { + res = r + if ( sum(abs(r)) !=0 ) + res = r/sum(abs(r)) + res +})) +colnames(signal_used) <- colnames(futures_data) +# Now create the returns of the strategy for each futures time series +strategy_by_future <- scrub(shift(signal_used,2)*futures_data) # use the signal from 2 days ago +# finally, this is our futures trend following strategy +trading_strategy = apply(strategy_by_future,1,sum) +names(trading_strategy) <- rownames(futures_data) +``` + + +### Reporting the performance results + +Let's see how this strategy does: +
+
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis',fig.align='center', fig.height=5,fig.width= 8, fig=TRUE} +pnl_plot(trading_strategy) +``` + +
+
+ +Here is how this strategy has performed during this period. +
+
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(cbind(rownames(pnl_matrix(trading_strategy)), round(pnl_matrix(trading_strategy),2))) +m1<-gvisTable(show_data,options=list(width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE)) +print(m1,'chart') +``` + +
+
+ +How does this compare with **existing CTA products** such as [this one from Societe Generale?](https://cib.societegenerale.com/fileadmin/indices_feeds/SG_CTA_Monthly_Report.pdf) (Note: one can easily achieve a correlation of more than 0.8 with this specific product - as well as with many other ones) + +![Compare our strategy with this product](societegenerale.png) + +
+ +#### Questions + +1. Can you describe in more detail what the code above does? +2. What happens if you use different moving average technical indicators in the code above? Please explore and report below the returns of a trading strategy you build. (Hint: check that the command line `MAfreq<-seq(10,250,by=20)` above does for example - but not only of course, the possibilities are endless) + +
+ +**Your Answers here:** + +1. The code generates a new variable (signal_used) which has the same format of futures_data. Then, it fills this variable with the "votes" of the moving averages of another new variable (MAfreq). These votes indicate the investment for each of the futures. + +2. If we change the values of the moving average (for example, giving a broader range), the returns will be slightly similar (although the trend is similar). + +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +signal_used = 0*futures_data +MAfreq<-seq(10,400,by=10) +for (iter in 1:length(MAfreq)) + signal_used = signal_used + sign(apply(futures_data,2, function(r) ma(r,MAfreq[iter]))) +signal_used = t(apply(signal_used,1,function(r) { + res = r + if ( sum(abs(r)) !=0 ) + res = r/sum(abs(r)) + res +})) +colnames(signal_used) <- colnames(futures_data) +strategy_by_future <- scrub(shift(signal_used,2)*futures_data) +trading_strategy = apply(strategy_by_future,1,sum) +names(trading_strategy) <- rownames(futures_data) +pnl_plot(trading_strategy) +``` + +
+ + +
+ +### A class competition + +Now you have seen how to develop some trading strategies that hedge funds have been using for centuries. Clearly this is only the very first step - as many of the online resources on technical indicators also suggest. Can you now explore more such strategies? How good a **futures trend following hedge fund strategy** can you develop? Let's call this.... a **class competition**! Explore as much as you can and report your best strategy as we move along the course... + +Here is for example something that can be achieved relatively easily... +
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis',fig.align='center', fig.height=5,fig.width= 8, fig=TRUE} +load("data/sample_strategy.Rdata") +pnl_plot(sample_strategy) +``` + +
+ +Here is how this strategy has performed during this period. +
+
+ +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data = data.frame(cbind(rownames(pnl_matrix(sample_strategy)), round(pnl_matrix(sample_strategy),2))) +m1<-gvisTable(show_data,options=list(width=1220, height=min(400,27*(nrow(show_data)+1)),allowHTML=TRUE)) +print(m1,'chart') +``` + +
+
+ +**Finally**: One can develop (shiny based) interactive versions of this report and deploy them using `shinyapps::deployApp('ExerciseSet2.Rmd')` (you need a [shinyapps.io](https://www.shinyapps.io) account for this). This is for example an [interactive version of this exercise.](https://inseaddataanalytics.shinyapps.io/ExerciseSet2/) + +
+
+ +As always, **have fun** + + + + + diff --git a/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.html b/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.html new file mode 100644 index 00000000..0daee929 --- /dev/null +++ b/Exercises/Exerciseset2/CopyOfExerciseSet2_Ana.html @@ -0,0 +1,9057 @@ + + + + + + + + + + + + + +Exercise Set 2: A $300 Billion Strategy + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +


+

The purpose of this exercise is to become familiar with:

+
    +
  1. Some time series analysis tools;
  2. +
  3. Correlation matrices and principal component analysis (PCA) (see readings of sessions 3-4);
  4. +
  5. More data manipulation and reporting tools (including Google Charts).
  6. +
+

As always, while doing this exercise we will also see how to generate replicable and customizable reports. For this purpose the exercise uses the R Markdown capabilities (see Markdown Cheat Sheet or a basic introduction to R Markdown). These capabilities allow us to create dynamic reports. For example today’s date is 2016-01-20 (you need to see the .Rmd to understand that this is not a static typed-in date but it changes every time you compile the .Rmd - if the date changed of course).

+

Before starting, make sure you have pulled the exercise set 2 souce code files on your github repository (if you pull the course github repository you also get the exercise set files automatically). Moreover, make sure you are in the directory of this exercise. Directory paths may be complicated, and sometimes a frustrating source of problems, so it is recommended that you use these R commands to find out your current working directory and, if needed, set it where you have the main files for the specific exercise/project (there are other ways, but for now just be aware of this path issue). For example, assuming we are now in the “Data Analytics R version/INSEADAnalytics” directory, we can do these:

+
getwd()
+setwd("Exercises/Exerciseset2/")
+list.files()
+

Note: as always, you can use the help command in Rstudio to find out about any R function (e.g. type help(list.files) to learn what the R function list.files does).

+

Let’s now see the exercise.

+

IMPORTANT: You should answer all questions by simply adding your code/answers in this document through editing the file ExerciseSet2.Rmd and then clicking on the “Knit HTML” button in RStudio. Once done, please post your .Rmd and html files in your github repository.

+
+
+

The Exercise: Introduction

+

For this exercise we will use the Futures’ daily returns to develop what is considered to be a “classic” hedge fund trading strategy, a futures trend following strategy. There is a lot written about this, so it is worth doing some online search about “futures trend following”, or “Managed Futures”, or “Commodity Trading Advisors (CTA)”. There is about $300 billion invested on this strategy today, and is considered to be one of the oldest hedge fund strategies. Some example links are:

+ +

Of course there are also many starting points for developing such a strategy (for example this R bloggers one (also on github), or the turtle traders website which has many resources.

+

In this exercise we will develop our own strategy from scratch.

+

Note (given today’s market conditions): Prices of commodities, like oil or gold, can be excellent indicators of the health of the economy and of various industries, as we will also see below.

+
+
+

Getting the Futures Data

+

There are many ways to get futures data. For example, one can use the Quandl package, or the turtle traders resources, or (for INSEAD only) get data from the INSEAD library finance data resources website. One has to pay attention on how to create continuous time series from underlying contracts with varying deliveries (e.g. see here ). Using a combination of the resources above, we will use data for a number of commodities.

+
+
+

Data description

+

Let’s load the data and see what we have.

+
source("helpersSet2.R")
+library(googleVis)
+load("data/FuturesTrendFollowingData.Rdata")
+


We have data from 2001-01-02 to 2015-09-24 of daily returns for the following 64 futures:

+


+
show_data = data.frame(colnames(futures_data))
+m1 <- gvisTable(show_data, options = list(showRowNumber = TRUE, width = 1920, 
+    height = min(400, 27 * (nrow(show_data) + 1)), allowHTML = TRUE, page = "disable"))
+print(m1, "chart")
+ + + + + + + +
+ +
+


+
+
+

Basic data analysis

+

Let’s see how these are correlated. Let’s also make it look nicer (than, say, what we did in Exercise Set 1), using Google Charts (see examples online, e.g. examples and the R package used used ).The correlation matrix is as follows (note that the table is “dynamic”: for example you can sort it based on each column by clicking on the column’s header)

+


+ + + + + + + +
+ +
+


+

We see quite high correlations among some of the futures. Does it make sense? Why? Do you see some negative correlations? Do those make sense?

+

Given such high correlations, we can try to see whether there are some “principal components” (see reading on dimensionality reduction). This analysis can also indicate whether all futures (the global economy!) are driven by some common “factors” (let’s call them “risk factors”).

+


+
Variance_Explained_Table_results <- PCA(futures_data, graph = FALSE)
+Variance_Explained_Table <- cbind(paste("component", 1:ncol(futures_data), sep = " "), 
+    Variance_Explained_Table_results$eig)
+Variance_Explained_Table <- as.data.frame(Variance_Explained_Table)
+colnames(Variance_Explained_Table) <- c("Component", "Eigenvalue", "Percentage_of_explained_variance", 
+    "Cumulative_percentage_of_explained_variance")
+ + + + + + + +
+ +
+


+

Here is the scree plot (see Sessions 3-4 readings):

+
eigenvalues <- Variance_Explained_Table[, 2]
+ + + + + + + +
+ +
+


+

Let’s now see how the 20 first (rotated) principal components look like. Let’s also use the rotated factors (note that these are not really the “principal component”, as explained in the reading on dimensionality reduction) and not show any numbers less than 0.3 in absolute value, to avoid cluttering. Note again that you can sort the table according to any column by clicking on the header of that column.

+
corused = cor(futures_data[, apply(futures_data != 0, 2, sum) > 10, drop = F])
+Rotated_Results <- principal(corused, nfactors = 20, rotate = "varimax", score = TRUE)
+Rotated_Factors <- round(Rotated_Results$loadings, 2)
+Rotated_Factors <- as.data.frame(unclass(Rotated_Factors))
+colnames(Rotated_Factors) <- paste("Component", 1:ncol(Rotated_Factors), sep = " ")
+
+sorted_rows <- sort(Rotated_Factors[, 1], decreasing = TRUE, index.return = TRUE)$ix
+Rotated_Factors <- Rotated_Factors[sorted_rows, ]
+Rotated_Factors[abs(Rotated_Factors) < 0.3] <- NA
+ + + + + + + +
+ +
+


+
+

Questions:

+
    +
  1. How many principal components (“factors”) do we need to explain at least 50% of the variance in this data?
  2. +
  3. What are the highest weights (in absolute value) of the first principal component portfolio above on the 64 futures?
  4. +
  5. Can we interpret the first 10 components? How would you call these factors?
  6. +
  7. Can you now generate the principal components and scree plot using only: a) the pre-crisis bull market years (e.g. only using the data between November 1, 2002, and October 1, 2007)? b) the financial crisis years (e.g. only using the data between October 1, 2007 and March 1, 2009), (Hint: you can select subsets of the data using for example the command `crisis_data = futures_data[as.Date(rownames(futures_data)) > “2007-10-01” & as.Date(rownames(futures_data)) < “2009-03-01”, ])
  8. +
  9. Based on your analysis in question 3, please discuss any differences you observe about the futures returns during bull and bear markets. What implications may these results have? What do the results imply about how assets are correlated during bear years compared to bull years?
  10. +
  11. (Extra - optional) Can you create an interactive (shiny based) tool so that we can study how the “risk factors” change ove time? (Hint: see Exercise set 1 and online resources on Shiny such as these Shiny lessons. Note however that you may need to pay attention to various details e.g. about how to include Google Charts in Shiny tools - so keep this extra exercise for later!).
  12. +
+


+

Your Answers here:

+
    +
  1. We need 6 factors (cumulative percentage = 52,29%)

  2. +
  3. In absolute values, the highest weights are 0,93 (5yr T-Notes US and 10yr T notes), 0,86 (treasury bonds), 0,85 (2 yr T-note US), 0,79, 0,74, …

  4. +
  5. If we see the first 10 factors, we can see that, since the overlapping of values between the factors are low, we can explain the variables with these 10 factors.

  6. +
  7. For the pre crisis data:

  8. +
+ + + + + +
+ +
+

For the crisis data:

+ + + + + + + +
+ +
+
    +
  1. During the bull (pre crisis) years we need more components to explain at least 50% of the data (again 6 factors, while during the crisis years we can explain it with 4 factors)
  2. +
+
+
+
+
+

A Simple Futures Trend Following Strategy

+

We can now develop a simple futures trend following trading strategy, as outlined in the papers in the Exercise Introduction above. There are about $300 billion invested in such strategies! Of course we cannot develop here a sophisticated product, but with some more work…

+

We will do the following:

+
    +
  1. Calculate a number of moving averages of different “window lengths” for each of the 64 futures - there are many so called technical indicators one can use. We will use the “moving average” function ma for this (try for example to see what this returns ma(1:10,2) ).
  2. +
  3. Add the signs (can also use the actual moving average values of course - try it!) of these moving averages (as if they “vote”), and then scale this sum across all futures so that the sum of their (of the sum across all futures!) absolute value across all futures is 1 (hence we invest $1 every day - you see why?).
  4. +
  5. Then invest every day in each of the 64 an amount that is defined by the weights calculated in step 2, using however the weights calculated using data until 2 days ago (why 2 days and not 1 day?) - see the use of the helper function shift for this.
  6. +
  7. Finally see the performance of this strategy.
  8. +
+

Here is the code.

+
signal_used = 0 * futures_data  # just initialize the trading signal to be 0
+# Take many moving Average (MA) Signals and let them 'vote' with their sign
+# (+-1, e.g. long or short vote, for each signal)
+MAfreq <- seq(20, 250, by = 30)
+for (iter in 1:length(MAfreq)) signal_used = signal_used + sign(apply(futures_data, 
+    2, function(r) ma(r, MAfreq[iter])))
+# Now make sure we invest $1 every day (so the sum of the absolute values of
+# the weights is 1 every day)
+signal_used = t(apply(signal_used, 1, function(r) {
+    res = r
+    if (sum(abs(r)) != 0) 
+        res = r/sum(abs(r))
+    res
+}))
+colnames(signal_used) <- colnames(futures_data)
+# Now create the returns of the strategy for each futures time series
+strategy_by_future <- scrub(shift(signal_used, 2) * futures_data)  # use the signal from 2 days ago
+# finally, this is our futures trend following strategy
+trading_strategy = apply(strategy_by_future, 1, sum)
+names(trading_strategy) <- rownames(futures_data)
+
+
+

Reporting the performance results

+

Let’s see how this strategy does:

+

+



+

Here is how this strategy has performed during this period.

+ + + + + + + +
+ +
+



+

How does this compare with existing CTA products such as this one from Societe Generale? (Note: one can easily achieve a correlation of more than 0.8 with this specific product - as well as with many other ones)

+

Compare our strategy with this product

+


+
+

Questions

+
    +
  1. Can you describe in more detail what the code above does?
  2. +
  3. What happens if you use different moving average technical indicators in the code above? Please explore and report below the returns of a trading strategy you build. (Hint: check that the command line MAfreq<-seq(10,250,by=20) above does for example - but not only of course, the possibilities are endless)
  4. +
+


+

Your Answers here:

+
    +
  1. The code generates a new variable (signal_used) which has the same format of futures_data. Then, it fills this variable with the “votes” of the moving averages of another new variable (MAfreq). These votes indicate the investment for each of the futures.

  2. +
  3. If we change the values of the moving average (for example, giving a broader range), the returns will be slightly similar (although the trend is similar).

  4. +
+
signal_used = 0 * futures_data
+MAfreq <- seq(10, 400, by = 10)
+for (iter in 1:length(MAfreq)) signal_used = signal_used + sign(apply(futures_data, 
+    2, function(r) ma(r, MAfreq[iter])))
+signal_used = t(apply(signal_used, 1, function(r) {
+    res = r
+    if (sum(abs(r)) != 0) 
+        res = r/sum(abs(r))
+    res
+}))
+colnames(signal_used) <- colnames(futures_data)
+strategy_by_future <- scrub(shift(signal_used, 2) * futures_data)
+trading_strategy = apply(strategy_by_future, 1, sum)
+names(trading_strategy) <- rownames(futures_data)
+pnl_plot(trading_strategy)
+

+


+
+
+
+
+

A class competition

+

Now you have seen how to develop some trading strategies that hedge funds have been using for centuries. Clearly this is only the very first step - as many of the online resources on technical indicators also suggest. Can you now explore more such strategies? How good a futures trend following hedge fund strategy can you develop? Let’s call this…. a class competition! Explore as much as you can and report your best strategy as we move along the course…

+

Here is for example something that can be achieved relatively easily…

+

+


+

Here is how this strategy has performed during this period.

+ + + + + + + +
+ +
+



+

Finally: One can develop (shiny based) interactive versions of this report and deploy them using shinyapps::deployApp('ExerciseSet2.Rmd') (you need a shinyapps.io account for this). This is for example an interactive version of this exercise.

+



+

As always, have fun

+
+ + +
+ + + + + + + +