[R] 複製－市盈率與後市關係之統計

最近好像少了Coding上的主題。上星期看到這篇文章：《以市盈率和股息率看後市統計》，正好想把這個作為R的練習。它統計在不同市盈率（P/E）之下，恆指在1個月以至1年之後是升還是跌。結論是：雖然恆指無論在高或低市盈率的短期(一個月)後，表現未見有特別；但當看一年或兩年後的話，則明顯升多跌少。而且市盈率越低，一年後及二年後上升的機會越大。

過程

想重複這個統計，需要的兩項資料分別是恆指收市價和恆指市盈率。在練習R時，收市價當然是用它的Library： "quantmod" 。背後的數據來源是Yahoo財經；因為今次需要的是月度資料（之後提到的恆指PE，要方便的話也是拿按月的資料），所以用compression="m"；再用"zoo"格式的as.yearmon(index( )) 紀錄這類月度的時間序列。

Prices = get.hist.quote(instrument="^hsi", start="1986-01-01",end="2016-05-31", 
  quote=c("Close","AdjClose"),provider="yahoo", 
  origin="1970-01-01",compression="m", retclass="zoo")

另一方面的市盈率資料，之前為了自己弄的Excel模型就有找過。可以用恆指公司網頁上的.xls檔案。因為在MacOS上沒有再找處理excel格式的套件，所以下載回來後，我先刪去了合併儲存格的大標題、改了日期格式、和另存成我熟悉的csv檔案。

這兩項數據有不同的日期表示方式，Yahoo用的是每月首個交易日，恆指公司用的卻是每月最後一日。所以先將兩項資料都用"zoo"格式的as.yearmon(index( )) 去轉換成配合這類月度資料的時間序列。將兩個數據合併，然後在所有數據中，再抽出所想要的時間片段。而與他的例子不同的是，我選擇了1987年1月至2016年5月的資料。

一開始只是計算某一個市盈率範圍在某個時間後的表現，這個可以比較簡單可以做到。然後就想要引伸到不同情況（5種市盈率範圍、5種觀察時期），還有要把每個結果整合到一個表格中來呈現，就只想到用nested-loop的方式去計算每一個數值，然後逐步組合成用於表達的Matrix。不過，其實這方法在R的世界中好像因為執行效率的考慮而儘量避免的，會比較推崇的是Vectorization，把數據一開始便弄成更多維度的矩陣去運算。

整個結果如下圖。

結果

首先是按市盈率看後市上市的機率，這個與他的結果差不多，都是反映市盈率這個指標比較適合有心看一年的「長線」投資者。以2016年5月的數字作例子解釋，當時恆指報：20815.09，市盈率：10.28倍。過去接近30年間，出現9-12倍PE的月份中，只有56%是一個月後上升的（98次當中有55次），但78%會在一年後是上升（92次當中有72次）。

下一樣想看的是：回報又是多少呢？這是當信號出現後，計算指定日子之後不論正／負的所有回報。繼續以上的例子，持貨1個月的回報只有+1.18%，但1年後平均有+14.28%。除了升市機會較高之外，更重要是的有更高回報的期待。

以上是R Consloe上的截圖，如果再整理出以下的表格應該可以更容易閱讀，並看到投資時間的影響。
比較 1 個月後：

	出現次數	上升次數	百份比	平均回報
PE<9	13	6	46%	0.20%
9<=PE<12	98	55	56%	1.18%
12<=PE<15	101	65	64%	1.34%
15<=PE<18	91	51	56%	0.97%
18<=PE<21	30	19	63%	0.97%
21<=PE	20	8	40%	-2.74%

比較 12 個月後：

	出現次數	上升次數	百份比	平均回報
PE<9	8	8	100%	56.00%
9<=PE<12	92	72	78%	14.00%
12<=PE<15	101	64	63%	15.00%
15<=PE<18	91	55	60%	6.43%
18<=PE<21	30	19	63%	-1.27%
21<=PE	20	8	40%	-12.25%

比較 24 個月後：

	出現次數	上升次數	百份比	平均回報
PE<9	8	8	100%	77.62%
9<=PE<12	80	66	83%	35.59%
12<=PE<15	101	74	73%	28.45%
15<=PE<18	91	54	59%	14.71%
18<=PE<21	30	20	67%	7.72%
21<=PE	20	3	15%	-17.96%

不過，最後還要看一看這張圖，再自行決定歷史的平均市盈率是否反映將來平均市盈率。據如林本利等的觀點，有人認為因為H股、價值重估、特殊收益等因素，現在的P/E與以往相比其實是被低估了，應該用去除這些因素後，一個數值較高的 「實際P/E」。所以，沿用一般定義下的P/E數據的話，會不會遇上近年的P/E 常態可能已經不同了的問題？

另外，市盈率代表的是市場上的股價和企業盈利的比率，如果將來公司的盈利下跌，其實可以在市盈率不變的情況中股價同時下跌。那麼，我們對未來的經濟展望如何？

R Source:

setwd("<MY PATH>")
library(quantmod);library(zoo);library(tseries);

#% Get Price from Yahoo
Prices = get.hist.quote(instrument="^hsi", start="1986-01-01",end="2016-05-31", quote=c("Close","AdjClose"),provider="yahoo", origin="1970-01-01",compression="m", retclass="zoo")

#% Get PE from HSI File
PE = read.zoo(read.table("hsi.csv", skip=0, sep=",", header=TRUE, na.strings="--"))

#% Combine information & choose window
index(Prices)  = as.yearmon(index(Prices))
index(PE)  = as.yearmon(index(PE))
myData = merge(Prices, PE[,1])
names(myData) = c("Close", "AdjClose", "PE")
pxData= window(myData, start = as.yearmon(as.Date("1986-01-01")), end = as.yearmon(as.Date("2016-05-31")))
 
#% Select
pxData$Signal = TRUE
pxData$Signal1 = pxData$PE<9 ignal2="pxData" pxdata="">=9 & pxData $PE<12 ignal3="pxData" pxdata="">=12 & pxData $PE<15 ignal4="pxData" pxdata="">=15 & pxData $PE<18 ignal5="pxData" pxdata="">=18 & pxData $PE<21 ignal6="pxData" pxdata="">=21

#% Process
Mat_cSignal = NULL; Mat_cSignalPos = NULL; Mat_rSignal = NULL; 
pxData$Closelag = pxData$Close

for (s in c("Signal1", "Signal2", "Signal3", "Signal4", "Signal5", "Signal6")){
  pxData$Signal = pxData[,s]
  #print(s); print("-------------------------")

  lagVec_cSignal = NULL; lagVec_cSignalPos = NULL; lagVec_rSignal = NULL;

  for (i in c(1,3,6,12,24)) {
    pxData$Closelag = NULL
    pxData$Closelag = lag(pxData$Close, i)
    pxData$Return=pxData$Closelag/pxData$Close - 1

    cSignal = length( which(pxData$Signal==1 & !is.na(pxData$Return)) )
    cSignalPos = length( which(pxData$Signal==1 & pxData$Return>0) )
    rSignal = mean(pxData[pxData$Signal==1,]$Return, na.rm=TRUE)

    #print(i); print(cSignal); print(cSignalPos); print(cSignalPos/cSignal); print(rSignal);

    lagVec_cSignal = c(lagVec_cSignal, cSignal)
    lagVec_cSignalPos = c(lagVec_cSignalPos, cSignalPos)
    lagVec_rSignal = c(lagVec_rSignal, rSignal)
  }
  Mat_cSignal = cbind(Mat_cSignal, lagVec_cSignal)
  Mat_cSignalPos = cbind(Mat_cSignalPos, lagVec_cSignalPos)
  Mat_rSignal = cbind(Mat_rSignal, lagVec_rSignal)

  #print("-------------------------")
}

#% Presentation & Demo
colNames = c("PE<9, "9<=PE<12", "12<=PE<15", "15<=PE<18", "18<=PE<21", "21<=PE")
rowNames = c("For 1 Month:", "For 3 Month:", "For 6 Month:", "For 12 Month:", "For 24 Month:")
colnames(Mat_cSignal) = colNames; colnames(Mat_cSignalPos) = colNames; colnames(Mat_rSignal) = colNames
rownames(Mat_cSignal) = rowNames; rownames(Mat_cSignalPos) = rowNames; rownames(Mat_rSignal) = rowNames
tail(pxData[,c(1,3)])
print("Number of Signal, since 1987:"); print(Mat_cSignal); 
print("Percentage of Rising Cases, after signal (%):"); print(Mat_cSignalPos/Mat_cSignal*100); 
print("Avg.Return of all Cases, after signal (%):"); print(Mat_rSignal*100); 
plot(pxData$PE, col='blue', lwd=2, ylab="P/E (HSI)" , xlab="Jan1987 to May2016")