diff --git a/CourseSessions/Sessions23/Session2inclass_BR.Rmd b/CourseSessions/Sessions23/Session2inclass_BR.Rmd new file mode 100644 index 00000000..dfb3c0b9 --- /dev/null +++ b/CourseSessions/Sessions23/Session2inclass_BR.Rmd @@ -0,0 +1,433 @@ +--- +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 responses of questions that are similiar ie. questions pertaining to interest in status symbols/brands or pertain to similiar habits ie. buying the greatest and latest boats and accessories or or questions relating to Boating as a passion and interest +
+2. These correlations imply that those who care about one thing care strongly about another thing that they are strongly correlated to. +
+ +
+ +### 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 helps us determine how many factors we should have. If the eigen value is greater than 1, then we should consider the component a factor. It finds the number of factors that are enough to be able to explain a signficant portion of the explained variance. The plot is a visual way of showing where the eigen values change from greater than 1 to less than 1. +
+2. The plot decreases because as the number of factors increases, the component's eigen values decreases until it nearly reaches one. The plotted line cannot increase. +
+3. We would prefer to see a smaller number of distinct factors as this will assist us in better defining clusers in the future. +
+ +#### 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 couple of factors relate to the most important factors as it relates to consumer attitudes and purchasing behaviors. Yes, this makes business sense because you can use the results to drive marketing decisions. +
+2. We should choose 5 factors based on the fact that 1) the eigen values are all greater than 1, 2) each of the factors represents a distinct characteristic, and 3) they are all important in driving marketing decisions. +
+3. Factor 1 refers to caring about the are look of the boat and user of boat as a status symbol. Factor 2 refers to adventure and loving the sense of freedom and activity that boating provides. Factor 3 refers to importance of brand and technology. Factor 4 refers to specialized knowledge about boats. Factor 5 refers to price sensitivity. +
+4. Data science goes hand in hand with business concepts (ie. marketing). You can't just look at numbers, but you have to understand business in order to interpret them and understand their relation to decision-making. +
+
+
+ +### 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:** +
1. There are 3 segments in the market. I chose 3 segments because 1) the dendrogram noted that there are three very distinct segments that we could classify; 2) when interpreting each segment, there appears very distinct customer characteristics that don't overlap; 3) there is a clear customer target that has emerged from this segmentation exercise. +
2. Segment 1 describes the customer that boats for social purposes and is not a boating enthusiast. Segment 2 describes the customer that refers to the customer who doesn't care about purchasing a boat now and Segment 3 relates to the customer who cares about boats and about being knowledgeable about boats. +
3.If I change the number of factors, I could create smaller clusters and do a microsegmentation exercise. However, those clusters may not be as distinct. +
4. The best customer target would be Segment 3. These customers care about a powerful boat that is affordable. They also tend to boat by themselves and to fix boats by themselves and are great resources about boating knowledge to their peers. They know so much about boating that they can identify what makes a great boat, but understands that high quality does not necessarily simply correlate just with price. +
+ +**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/Exercises/Exerciseset1/DataSet1.Rdata b/Exercises/Exerciseset1/DataSet1.Rdata new file mode 100644 index 00000000..ba4222cd 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..22a14f45 100644 --- a/Exercises/Exerciseset1/ExerciseSet1.Rmd +++ b/Exercises/Exerciseset1/ExerciseSet1.Rmd @@ -41,6 +41,7 @@ We download daily prices (open, high, low, close, and adjusted close) and volume ```{r eval = TRUE, echo=TRUE, error = FALSE, warning=FALSE,message=FALSE,results='asis'} source("helpersSet1.R") source("dataSet1.R") +source("dataSet2.R") ``` For more information on downloading finance data from the internet as well as on finance related R tools see these starting points (there is a lot more of course available): @@ -55,12 +56,12 @@ For more information on downloading finance data from the internet as well as on #### Optional Question -1. Can you find some interesting finance related R package or github repository? -**Your Answers here:** -
-
- +1. Can you find some interesting finance related R package or github repository?
+1. *Several R finance package are interesting including:*
+*- Rblpapi: a R Interface to Bloomberg*
+*- FinCal: a package to calculate the time Value of Money, the time Series Analysis and computational Finance*
+*- jrvFinance: a basic finance tool to calculate NPV/IRR/Annuities/Bond-Pricing or Black Scholes*
### Part I: Statistics of S&P Daily Returns @@ -86,11 +87,16 @@ pnl_plot(SPY) **Your Answers here:**
-
-
-
+ +*line 96: *ExerciseSet1.Rmd:write.csv(StockReturns[1:20,c("SPY","**AAPL**")], file = "twentydays.csv", row.names = TRUE, col.names = TRUE)
+*line 8: *dataSet1.R:mytickers = c("SPY", "**AAPL**") # Other tickers for example are "GOOG", "GS", "TSLA", "FB", "MSFT",
+ +
+*Cumulative returns:* `r round(100*sum(StockReturns[,2]),1)`%
+*Average daily returns:* `r round(100*mean(StockReturns[,2]),3)`%
+*Standard deviation of the daily returns of Apple:* `r round(100*sd(StockReturns[,2]),3)`%
### Part II: Simple Matrix Manipulations @@ -125,13 +131,15 @@ We can also transpose the matrix of returns to create a new "horizontal" matrix. 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 can use `nrow(transposedData)` and `ncol(transposedData)`. +2. I assume that we are asked to re-draw the plot using the transposed matrix: +```{r eval=FALSE, echo=FALSE, comment=NA, warning=FALSE, message=FALSE,results='asis',fig.align='center', fig=TRUE} +portfolio2 = apply(transposedData,2,mean) +names(portfolio2) <- colnames(transposedData) +pnl_plot(portfolio2) +``` +
### Part III: Reproducibility and Customization @@ -142,14 +150,8 @@ This is an important step and will get you to think about the overall process on 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:** -
-
-
-
-
+All we need to change is the value allocated to the variable `startDate = "2001-01-01"` in Line 11.
### Part IV: Read/Write .CSV files @@ -177,13 +179,12 @@ Try it! myData + StockReturns[1:40,]) ``` -**Your Answers here:** -
-
-
-
-
+ +1. We get as output a number (20). Each row in myData is different from the corresponding row in StockReturns (hence the test is true and yields a result of 1). This makes senses because myData contains one row with a string whereas StockReturns containts three rows. +However, if we use as separator in the read.csv function a comma instead of a semicolon, we end up with matrices of the same dimensions, but there is an offset of one row for some reason. THen we get a result of 60 (20x3). If we were able to resolve the offset issue, we'd get a result of zero, i.e. all the cells would match. +2. We will get an error because we are summing matrices with different types of arguments. +
### Extra Question diff --git a/Exercises/Exerciseset1/ExerciseSet1.html b/Exercises/Exerciseset1/ExerciseSet1.html new file mode 100644 index 00000000..e0af4dd7 --- /dev/null +++ b/Exercises/Exerciseset1/ExerciseSet1.html @@ -0,0 +1,254 @@ + + + + + + + + + + + + + +Exercise Set 1 + + + + + + + + + + + + + + + + + + + + + +
+ + + + + +


+

The purpose of this exercise is to become familiar with:

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

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-26 (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). For example, assuming we are now in the “MYDIRECTORY/INSEADAnalytics” directory, we can do these:

+
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")
+

[1] “ticker SPY …” [1] “ticker AAPL …”

+
source("dataSet2.R")
+

[1] “ticker SPY …” [1] “ticker AAPL …”

+

For more information on downloading finance data from the internet as well as on finance related R tools see these starting points (there is a lot more of course available):

+ +
+

Optional Question

+
    +
  1. Can you find some interesting finance related R package or github repository? +
    2. +
  2. Several R finance package are interesting including:
    - Rblpapi: a R Interface to Bloomberg
    - FinCal: a package to calculate the time Value of Money, the time Series Analysis and computational Finance
    - jrvFinance: a basic finance tool to calculate NPV/IRR/Annuities/Bond-Pricing or Black Scholes
    +
    3. +
+
+
+
+

Part I: Statistics of S&P Daily Returns

+

We have 2783 days of data, starting from 2005-01-04 until 2016-01-25. Here are some basic statistics about the S&P returns:

+
    +
  1. The cumulative returns of the S&P index during this period is 88.7%.
    2. +
  2. The average daily returns of the S&P index during this period is 0.032%;
    3. +
  3. The standard deviation of the daily returns of the S&P index during this period is 1.259%;
    4. +
+

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. +
  2. What are the cumulative, average daily returns, and the standard deviation of the daily returns of Apple in the same period?
    3. +
  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?
    4. +
+

Your Answers here:

+
+

line 96: ExerciseSet1.Rmd:write.csv(StockReturns[1:20,c(“SPY”,“AAPL”)], file = “twentydays.csv”, row.names = TRUE, col.names = TRUE)
line 8: dataSet1.R:mytickers = c(“SPY”, “AAPL”) # Other tickers for example are “GOOG”, “GS”, “TSLA”, “FB”, “MSFT”,

+
+Cumulative returns: 383.2%
Average daily returns: 0.138%
Standard deviation of the daily returns of Apple: 2.19%
+
+
+
+
+

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. +
  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?
    3. +
+
+
    +
  1. We can use nrow(transposedData) and ncol(transposedData).
    2. +
  2. I assume that we are asked to re-draw the plot using the transposed matrix:
    3. +
+
+
+
+
+

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. +
  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? +
    +All we need to change is the value allocated to the variable startDate = "2001-01-01" in Line 11. +
    3. +
+
+
+
+

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. +
  2. (Extra exercise) What do you think will happen if you now run this command, and why:
    3. +
+
myData + StockReturns[1:40,])
+
+
    +
  1. We get as output a number (20). Each row in myData is different from the corresponding row in StockReturns (hence the test is true and yields a result of 1). This makes senses because myData contains one row with a string whereas StockReturns containts three rows. However, if we use as separator in the read.csv function a comma instead of a semicolon, we end up with matrices of the same dimensions, but there is an offset of one row for some reason. THen we get a result of 60 (20x3). If we were able to resolve the offset issue, we’d get a result of zero, i.e. all the cells would match.
    2. +
  2. We will get an error because we are summing matrices with different types of arguments.
    3. +
+
+
+
+
+

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])
+})
+


+
+
+
+
+

Endless explorations (optional homework)

+

This is a recent research article that won an award in 2016. Can you implement a simple strategy as in Figure 1 of this paper? You may find these R commands useful: names, which, str_sub,diff,as.vector, length, pmin, pmax, sapply, lapply,Reduce,unique, as.numeric, %in% A Simple Trading Startegy

+

What if you also include information about bonds? (e.g. download the returns of the the ETF with ticker “TLT”) Is there any relation between stocks and bonds?

+

Have fun

+
+ + +
+ + + + + + + + diff --git a/Exercises/Exerciseset1/dataSet2.R b/Exercises/Exerciseset1/dataSet2.R new file mode 100644 index 00000000..c1917439 --- /dev/null +++ b/Exercises/Exerciseset1/dataSet2.R @@ -0,0 +1,51 @@ + +# rm(list=ls()) # Clean up the memory, if we want to rerun from scratch +# source("helpersSet1.R") + +getdata.fromscratch = 1 + +website_used = "yahoo" # can be "yahoo" or other ( see help(getSymbols) ). Depending on the website we may need to change the stock tickers' representation +mytickers = c("SPY", "AAPL") # Other tickers for example are "GOOG", "GS", "TSLA", "FB", "MSFT", +startDate = "2005-01-01" + +if (getdata.fromscratch){ + # Get SPY first, to get all trading days + tmp<-as.matrix(try(getSymbols(Symbols="AAPL",from = startDate,src = website_used, auto.assign=FALSE))) + ApplStockPrices=matrix(rep(0,nrow(tmp)*length(mytickers)), ncol=length(mytickers)) + colnames(ApplStockPrices)<-mytickers; + rownames(ApplStockPrices)<-rownames(tmp) + ApplStockVolume=ApplStockPrices + ApplStockPrices[,1] <- tmp[,6] + + for (ticker_index in 1:length(mytickers)){ + ticker_to_get = mytickers[ticker_index] + print(paste("\nDownloading ticker ", ticker_to_get, " ...")) + tmpdata<-as.matrix(try(getSymbols(Symbols=ticker_to_get,from = startDate,auto.assign=FALSE))) + if (!inherits(tmpdata, "try-error")) + { + therownames=intersect(rownames(tmpdata),rownames(ApplStockPrices)) + tmpdata[is.na(tmpdata)] <- 0 + ApplStockPrices[therownames,ticker_index]<-tmpdata[therownames,6] # adjusted close price + ApplStockVolume[therownames,ticker_index]<-tmpdata[therownames,5] # shares volume for now - need to convert to dollars later + } else { + cat(ticker_to_get," NOT found") + } + } + # Get the daily returns now. Use the simple percentage difference approach + AaplStockReturns= ifelse(head(ApplStockPrices,-1)!=0, (tail(StockPrices,-1)-head(StockPrices,-1))/head(StockPrices,-1),0) # note that this removes the first day as we have no way to get the returns then! + rownames(AaplStockReturns)<-tail(rownames(ApplStockPrices),-1) # adjust the dates by 1 day now + + # Now remove the first day from the other data, too + ApplStockPrices = StockPrices[rownames(AaplStockReturns),] + ApplStockVolume = ApplStockPrices[rownames(AaplStockReturns),] + colnames(ApplStockPrices)<-mytickers + colnames(ApplStockVolume)<-mytickers + colnames(AaplStockReturns)<-mytickers + + save(AaplStockReturns,ApplStockPrices,ApplStockVolume, file = "DataSet1.Rdata") +} else { + load("DataSet1.Rdata") +} + + + diff --git a/Exercises/Exerciseset1/twentydays.csv b/Exercises/Exerciseset1/twentydays.csv new file mode 100644 index 00000000..860a28c9 --- /dev/null +++ b/Exercises/Exerciseset1/twentydays.csv @@ -0,0 +1,21 @@ +"","SPY","AAPL" +"2005-01-04",-0.0122194637127275,0.0102701088517504 +"2005-01-05",-0.00690061615719681,0.00875830945932705 +"2005-01-06",0.00508430540389297,0.000775098385044691 +"2005-01-07",-0.00143324918157237,0.0728119031851786 +"2005-01-10",0.00472810453179597,-0.00418786145491284 +"2005-01-11",-0.00689074776905752,-0.063805107443653 +"2005-01-12",0.00330004252922032,0.0139406693608039 +"2005-01-13",-0.00801211937901914,0.066299907249727 +"2005-01-14",0.0052711748835299,0.00573060038933524 +"2005-01-18",0.0104026005444693,0.00641032505910421 +"2005-01-19",-0.0104628791052045,-0.0108989459656001 +"2005-01-20",-0.0060903540542225,0.00829999885962422 +"2005-01-21",-0.00612767375693113,0.0004257211545643 +"2005-01-24",-0.00196947533623822,0.00383049984844807 +"2005-01-25",0.00283135035218729,0.0182305207326887 +"2005-01-26",0.00299458297525282,0.0027759262791443 +"2005-01-27",0.00170602341346994,0.00539788429906169 +"2005-01-28",0,0.0184471453404301 +"2005-01-31",0.00621649614031668,0.0394701515667315 +"2005-02-01",0.00634732427396169,0.00819242781131764 diff --git a/Exercises/Exerciseset2/ExerciseSet2_BR.Rmd b/Exercises/Exerciseset2/ExerciseSet2_BR.Rmd new file mode 100644 index 00000000..b67034d6 --- /dev/null +++ b/Exercises/Exerciseset2/ExerciseSet2_BR.Rmd @@ -0,0 +1,696 @@ +--- +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. 6 components (cumulative % variance explained is 52%) +
2. 5yr and 10yr Treasury notes have highest weights of 0.93 +
3. The first ten components are agricultural commodities +
4. a) pre-crisis bull market years (e.g. only using the data between November 1, 2002, and October 1, 2007) +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +precrisis_data = futures_data[as.Date(rownames(futures_data)) > "2002-11-01" & as.Date(rownames(futures_data)) < "2007-10-01", ] + +Variance_Explained_Table_results_precrisis <- PCA(precrisis_data, graph=FALSE) +Variance_Explained_Table_precrisis <- cbind(paste("component",1:ncol(precrisis_data),sep=" "),Variance_Explained_Table_results_precrisis$eig) +Variance_Explained_Table_precrisis <- as.data.frame(Variance_Explained_Table_precrisis) +colnames(Variance_Explained_Table_precrisis) <- c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") + +eigenvalues_precrisis <- Variance_Explained_Table_precrisis[,2] +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data_precrisis = data.frame(Variance_Explained_Table_precrisis) +m1_precrisis <-gvisTable(show_data_precrisis,options=list(width=1920, height=min(400,27*(nrow(show_data_precrisis)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1_precrisis,'chart') +``` + +```{r Fig2, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df_precrisis <- cbind(as.data.frame(eigenvalues_precrisis), c(1:length(eigenvalues_precrisis)), rep(1, length(eigenvalues_precrisis))) +colnames(df_precrisis) <- c("eigenvalues", "components", "abline") +Line_precrisis <- gvisLineChart(as.data.frame(df_precrisis), 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_precrisis'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line_precrisis, 'chart') +``` + +
b) the financial crisis years (e.g. only using the data between October 1, 2007 and March 1, 2009) +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +crisis_data = futures_data[as.Date(rownames(futures_data)) > "2007-10-01" & as.Date(rownames(futures_data)) < "2009-03-01", ] + +Variance_Explained_Table_results_crisis <- PCA(crisis_data, graph=FALSE) +Variance_Explained_Table_crisis <- cbind(paste("component",1:ncol(crisis_data),sep=" "),Variance_Explained_Table_results_crisis$eig) +Variance_Explained_Table_crisis <- as.data.frame(Variance_Explained_Table_crisis) +colnames(Variance_Explained_Table_crisis) <- c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") + +eigenvalues_crisis <- Variance_Explained_Table_crisis[,2] +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data_crisis = data.frame(Variance_Explained_Table_crisis) +m1_crisis <-gvisTable(show_data_crisis,options=list(width=1920, height=min(400,27*(nrow(show_data_crisis)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1_crisis,'chart') +``` + +```{r Fig3, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df_crisis <- cbind(as.data.frame(eigenvalues_crisis), c(1:length(eigenvalues_crisis)), rep(1, length(eigenvalues_crisis))) +colnames(df_crisis) <- c("eigenvalues", "components", "abline") +Line_crisis <- gvisLineChart(as.data.frame(df_crisis), 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_crisis'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line_crisis, 'chart') +``` + +
+
5. It takes less components to explain 50% of the variance during the crisis years, suggesting that during bear years the returns of agricultural trend with the overall market more than bull years. +
+ +
+ +### 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(10,250,by=20) +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:** +
+
+
+
+ + +
+ +### 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** + + + + + +--- +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. 6 components (cumulative % variance explained is 52%) +
2. 5yr and 10yr Treasury notes have highest weights of 0.93 +
3. The first ten components are agricultural commodities +
4. a) pre-crisis bull market years (e.g. only using the data between November 1, 2002, and October 1, 2007) +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +precrisis_data = futures_data[as.Date(rownames(futures_data)) > "2002-11-01" & as.Date(rownames(futures_data)) < "2007-10-01", ] + +Variance_Explained_Table_results_precrisis <- PCA(precrisis_data, graph=FALSE) +Variance_Explained_Table_precrisis <- cbind(paste("component",1:ncol(precrisis_data),sep=" "),Variance_Explained_Table_results_precrisis$eig) +Variance_Explained_Table_precrisis <- as.data.frame(Variance_Explained_Table_precrisis) +colnames(Variance_Explained_Table_precrisis) <- c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") + +eigenvalues_precrisis <- Variance_Explained_Table_precrisis[,2] +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data_precrisis = data.frame(Variance_Explained_Table_precrisis) +m1_precrisis <-gvisTable(show_data_precrisis,options=list(width=1920, height=min(400,27*(nrow(show_data_precrisis)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1_precrisis,'chart') +``` + +```{r Fig2, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df_precrisis <- cbind(as.data.frame(eigenvalues_precrisis), c(1:length(eigenvalues_precrisis)), rep(1, length(eigenvalues_precrisis))) +colnames(df_precrisis) <- c("eigenvalues", "components", "abline") +Line_precrisis <- gvisLineChart(as.data.frame(df_precrisis), 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_precrisis'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line_precrisis, 'chart') +``` + +
b) the financial crisis years (e.g. only using the data between October 1, 2007 and March 1, 2009) +```{r echo=TRUE, eval=TRUE, comment=NA, warning=FALSE,error=FALSE, message=FALSE, prompt=FALSE, tidy=TRUE} +crisis_data = futures_data[as.Date(rownames(futures_data)) > "2007-10-01" & as.Date(rownames(futures_data)) < "2009-03-01", ] + +Variance_Explained_Table_results_crisis <- PCA(crisis_data, graph=FALSE) +Variance_Explained_Table_crisis <- cbind(paste("component",1:ncol(crisis_data),sep=" "),Variance_Explained_Table_results_crisis$eig) +Variance_Explained_Table_crisis <- as.data.frame(Variance_Explained_Table_crisis) +colnames(Variance_Explained_Table_crisis) <- c("Component","Eigenvalue", "Percentage_of_explained_variance", "Cumulative_percentage_of_explained_variance") + +eigenvalues_crisis <- Variance_Explained_Table_crisis[,2] +``` + +```{r echo=FALSE, comment=NA, warning=FALSE, message=FALSE, results='asis'} +show_data_crisis = data.frame(Variance_Explained_Table_crisis) +m1_crisis <-gvisTable(show_data_crisis,options=list(width=1920, height=min(400,27*(nrow(show_data_crisis)+1)),allowHTML=TRUE,page='disable'),formats=list(Eigenvalue="#.##",Percentage_of_explained_variance="#.##",Cumulative_percentage_of_explained_variance="#.##")) +print(m1_crisis,'chart') +``` + +```{r Fig3, echo=FALSE, comment=NA, results='asis', message=FALSE, fig.align='center', fig=TRUE} +df_crisis <- cbind(as.data.frame(eigenvalues_crisis), c(1:length(eigenvalues_crisis)), rep(1, length(eigenvalues_crisis))) +colnames(df_crisis) <- c("eigenvalues", "components", "abline") +Line_crisis <- gvisLineChart(as.data.frame(df_crisis), 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_crisis'}]", series="[{color:'green',pointSize:3, targetAxisIndex: 0}]")) +print(Line_crisis, 'chart') +``` + +
+
5. It takes less components to explain 50% of the variance during the crisis years, suggesting that during bear years the returns of agricultural trend with the overall market more than bull years. +
+ +
+ +### 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(10,250,by=20) +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:** +
+
+
+
+ + +
+ +### 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/INSEADAnalytics.Rproj b/INSEADAnalytics.Rproj new file mode 100644 index 00000000..8e3c2ebc --- /dev/null +++ b/INSEADAnalytics.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX