作者：謝佳標

中國R語言大會講師，高階資料分析師，8年以上資料探勘建模工作實戰經驗

https://ask.hellobi.com/blog/xiejiabiao/4288

利用主成分分析構造你個人的股市指數，然後分析你的私家指數和該股市常用官方股票指數的相關性。

接用ML_for_Hackers-master 書中的資料。

> prices<-read.csv("stock_prices.csv")
> prices[1,]
        Date Stock Close
1 2011-05-25   DTE 51.12
> # 原始資料集並不是我們喜歡使用的格式，因此需要進行預處理。
> # 第一步，把資料集中的時間戳轉換為正確編碼的日期變數。這要用到lubridate包中的ymd函式
> # install.packages("lubridate")
> library(lubridate)
> prices<-transform(prices,Date=ymd(Date))
> # 一旦完成這一步，就能適用reshape函式庫中的cast函式
> library(reshape)

Attaching package: ‘reshape’

The following object is masked from ‘package:lubridate’:

    stamp

> date.stock.matrix<-cast(prices,Date~Stock,value="Close")
> which(complete.cases(date.stock.matrix)==F) # 22 875條記錄有缺失值
[1]  22 875
> # 分析了這個生成結--巨大的日期-股票矩陣之後，我們注意到缺失了一些元素。
> date.stock.matrix[22,];date.stock.matrix[875,]
         Date ADC AFL ARKR AZPN CLFD DDR DTE ENDP FLWS FR GMXR GPC HE ISSC ISSI KSS MTSC
22 2002-02-01  NA  NA   NA   NA   NA  19  NA   NA   NA NA   NA  NA NA   NA   NA  NA   NA
   NWN ODFL PARL RELV SIGM STT TRIB UTR
22  NA   NA   NA   NA   NA  NA   NA  NA
          Date  ADC   AFL  ARKR AZPN CLFD DDR   DTE  ENDP FLWS    FR  GMXR   GPC    HE
875 2005-06-22 30.4 43.49 26.56 5.76 1.47  NA 46.89 25.88 7.23 41.45 13.45 42.76 27.21
     ISSC ISSI   KSS  MTSC   NWN  ODFL  PARL RELV SIGM   STT TRIB   UTR
875 35.42 7.22 56.06 34.54 36.87 27.74 28.96 10.3 8.12 49.22 6.69 49.98
> # 因此回到最初的prices資料集，刪除那些缺失元素的資料，再執行cast函式：
> prices<-subset(prices,Date!=ymd('2002-02-01'))
> prices<-subset(prices,Stock !='DDR')
> date.stock.matrix<-cast(prices,Date~Stock,value="Close")
> which(complete.cases(date.stock.matrix)==F)  
#0
integer(0)
> # 接下來可以適用cor函式來找到這個矩陣中所有數字列之間的相關性。然後將相關性矩陣轉換成一個數值向量，並且畫一個相關性密度圖，
> # 以此來獲得兩個直觀認識：a）相關性的均值；b）低相關性出現的頻率。
> cor.matrix<-cor(date.stock.matrix[,2:ncol(date.stock.matrix)])
> correlations<-as.numeric(cor.matrix)
> library(ggplot2)
> ggplot(data.frame(Correlation=correlations),
+        aes(x=Correlation,fill=1))+
+   geom_density()+opts(legend.position='none')
 

> #正如密度圖所示，大部分相關性是正數，因此PCA適合用於這份資料集.
> # 我們適用princomp函式來執行PCA：
> pca<-princomp(date.stock.matrix[,2:ncol(date.stock.matrix)])
> # 我們只對第一主成份感興趣，所以只把pca載荷的第一列提取出來：
> principal.component<-pca$loadings[,1]
> # 完成這些之後，我們可以分析載荷的密度圖，直觀地瞭解第一主成份是如何形成的。
> loadings<-as.numeric(principal.component)
> ggplot(data.frame(Loadings=loadings),
+        aes(x=Loadings,fill=1))+
+   geom_density()+opts(legend.position="none")
 

> # 這個結果有點讓人疑惑，因為載荷有一個相當不錯的分佈，但是幾乎全是負數。它實際上是個很小的麻煩，我們用一行程式碼就能解決。
> # 到目前為止我們獲得了主成分，接下來可以把這些資料總結成一列了。可以使用predict函式完成這個目標：
> market.index<-predict(pca)[,1]
> # 如何才能知道這些預測值的效果呢？幸運的是，對這個例項我們可以很容易地判斷結果好壞

> # ，因為可以把結果和著名的市場指數做比較。在本章中，我們用道瓊斯指數（Down Jones Index DJI).
> dji.prices<-read.csv("DJI.csv")
> dji.prices<-transform(dji.prices,Date=ymd(Date))
> # 因為使用整個DJI執行的時間比我們預想的要長很多，所以需要取一個它的自己，僅僅獲得我們感興趣的那些日期。
> dji.prices<-subset(dji.prices,Date>ymd('2001-12-31'))
> dji.prices<-subset(dji.prices,Date !=ymd('2002-02-01'))
> # 然後，提取DJI中我們感興趣的部分，也就是每日收盤價格和我們記錄過的那些日期。
> # 因為它們的順序和我們現在的資料集相反，用rev函式反轉它們即可：
> dji<-with(dji.prices,rev(Close))
> dates<-with(dji.prices,rev(Date))
> # 現在我們可以繪製一些簡單的圖，將使用PCA生成的市場指數和DJI相比較：
> comparison<-data.frame(Date=dates,MarketIndex=market.index,DJI=dji)
> ggplot(comparison,aes(x=MarketIndex,y=DJI))+
+   geom_point()+geom_smooth(method="lm",se=FALSE)
 

> #從圖可以看出，那些之前看上去煩人的負載荷，真的成為了麻煩的源頭：我們的指數和DJI付相關。
> # 但是，我們可以很容易地解決這個麻煩。只需要對指數乘以-1，即可生成一個和DJI正相關的指數
> comparison<-transform(comparison,MarketIndex=-1*MarketIndex)
> # 現在可以再嘗試一次進行比較：
> ggplot(comparison,aes(x=MarketIndex,y=DJI))+
+   geom_point()+geom_smooth(method="lm",se=FALSE)
+   geom_point()+geom_line(
)

> #如圖，我們已經修正了指數的方向，並且它看上去和DJI真的很匹配。
> # 剩下的最後一件事情，就是獲得我們的指數隨著時間推移與DJI的趨勢保持一直的程度。
> # 首先，使用melt函式獲得一個數據框，它可以很容易地一次性對兩個指標進行視覺化。
> # 然後，我們對每個指數活出一條以日期為x軸，以價格為y軸的線。
> alt.comparison<-melt(comparison,id.vars="Date")
> names(alt.comparison)<-c("Date","Index","Price")
> ggplot(alt.comparison,aes(x=Date,y=Price,group=Index,color=Index))+
+   geom_point()+geom_line()

> # 這一次結果並不是很好，因為DJI都是很高的值，而我們的指數都是很小的值，但是可以使用scale函式解決這個問題。
> comparison$MarketIndex<-scale(comparison$MarketIndex)
> comparison$DJI<-scale(comparison$DJI)
> alt.comparison<-melt(comparison,id.vars="Date")
> names(alt.comparison)<-c("Date","Index","Price")
> ggplot(alt.comparison,aes(x=Date,y=Price,group=Index,color=Index))+
+   geom_point()+geom_line()

> # 看上去與DJI的趨勢保持得相當好。總之，用PCA真的能夠產生一副股票價格的趨勢圖。