時間數列 (二) – R 實作

R提供一個特殊的資料型態 ts, 作為Time Series基本類別，並提供相關方法, 作法與 matrix 類似, 這真是太厲害了，完全是OOP的標準作法，讓我對R的發明者更欽佩三分。

以下就是一個簡的例子
> library(MASS) # 匯入測試資料
> plot(huron,type="l") # 製作線圖
> setwd("j:/r/my") # 設定工作目錄
在繪圖視窗按滑鼠右鍵, 選【儲存BMP】，再貼到【小畫家】中存檔
> y2=huron$mean.height*2 # 產生應變數數列
> y=ts(y2) # 將data frame資料轉為time series資料型態
> x=huron$mean.height # 產生自變數數列
> x=ts(x)
> x2=lag(x,-1) # 產生第二個自變數數列
> d=ts.union(y,x,x2) # 組成一個ts資料來源
> lm(y~x+x2, data=d) #產生線性模型
lm(formula = y ~ x + x2, data = d)
Coefficients:
(Intercept) x x2
2.572e-13 2.000e+00 -4.273e-17

時間數列 (一) – 基本理論

一般迴歸分析(Regression)是假設資料誤差是完全隨機的，但現實社會中的現象，卻往往不是如此，例如股價，它的波動會受到前一日或更前幾天的影響，因此，要建構這類的模型，我們就可以應用時間數列(Time Series)，來分析或預測未來的走勢。

在談到時間數列前，我們最好先了解迴歸模型，以簡單線性迴歸為例:

y=a+bx+ε

其中ε為誤差, 假設他服從常態分數，平均數為0，變異數為σ²，且ε_i及ε_j無關聯

之後我們用誤差的最小平方法(Least Square)，求出a及b，進而得到

1. 總平方和(SST) = 殘差平方和(SSE) + 迴歸平方和(SSR)
   
2. R² = SSR / SST，R²是檢定x與y的相關係數
   
3. 檢定一元迴歸可用自由度=n-2的t分配，多元迴歸則用F分配
   

   F = (SSR / (k -1)) / (SSE / (n – k))
   

   其中k-1為自變數個數
   

時間數列與迴歸分析不同點就是假設誤差不是完全隨機的, 而是自我相關的，即前期的誤差會影響後期的誤差

ε_t = ρε_{t -1} + a_t

其中ρ絕對值<1，a_t稱為干擾項(Disturbance Term)，服從常態分數，平均數為0，變異數為σ²

1. 檢定ρ是否為0，一階(誤差只與前一期有關)可用Dubin-Watson統計量檢定，多階可用ARIMA統計量。
   
2. ARIMA： Auto Regression Integrated Moving Average，為Box/Jenkins提出
   

此外，若考慮不是σ²固定常數，也就是異質變異數(Heteroscendasticity)，可用加權最小平方法處理，這就是ARCH或GARCH模型。

註:此文主要是節錄自林茂文先生所注之時間數列之分析與預測一書。

迴歸(Regression)

重拾統計, 真是令人緬懷當年學生生活.

就從迴歸作起吧.

1. 載入測試資料, 也可以自己輸入
   > load("j:/r/programs/testdata/usingr.rdata")
2. 使用 moths 變數, 先輸入moths, 查看資料
3. 以moths$meters 為應變數, moths$A、moths$P 為自變數
   > lm(moths$meters ~ moths$A+moths$P)
4. 得到簡單迴歸的截距(Interception)及斜率(Slope)
   Coefficients:(Intercept) moths$A moths$P
   45.1017 0.3712 -0.6005

R 的編輯器

R提供較簡易的介面, 方便使用者輸入資料, 試試看下列方式吧.

1. b=de(a) # de 為 data.entry的簡稱, b 資料型態為 list
2. b=fix(a) # fix 為指令編輯視窗
3. b=edit(a) # 可輸入數學式子, 也可呼叫外部編輯器

執行外部程式

R可以下指令執行外部程式, 執行DOS指令時須加command /c, 用法如2

1. system("c:/windows/notepad.exe")
2. system(paste(Sys.getenv("COMSPEC"),"/c", "dir"))

相反的, 外部程式要執行R的程式也有兩種

1. R xxx.r # xxx.r 為R的程式碼
2. 使用com介面, 須加裝 RSrv250_pl1.exe, 內附範例可參考

R 變數類型

R的變數類型有下列幾種:

1. 向量: 為預設類型, 即一維陣列
2. Data frame: 具有欄位名稱的二維陣列
3. List: 可儲存不同類型的物件集合, 例如summary(), 要轉成向量, 可用unlist()
4. S3: 可具有class的list, 如果list被函數呼叫時, R會使用class的方法取代標準的方法, 請看下例:
   # 取三個常態分配隨機亂數
   h=list(a=rnorm(3), b= "it cannot be printed ")
   # 設定h變數的class名稱為MyClass
   class(h)= "MyClass "
   # 實作MyClass函數
   print.MyClass=function(x){cat("A = ", x$a, "n ")}
   # 列印h, 發現R以MyClass取代print
   print(h)
5. S4: 類似S3, 有點複雜, 以後再說明
6. 日期變數: 要使用日期函數處理, 必須先把變數轉為日期
   l 如為字串變數且格式為yyyy-mm-dd, 可用as.Date轉換
   l 如為數字變數且格式為yyyymmdd, 須先用as.character轉為文字再使用as.Date轉換
   > a=20081028
   > as.Date(as.character(a), format="%Y%m%d")
   [1] "2008-10-28"
   > as.Date(as.character(a), format="%Y%m%d")+1
   [1] "2008-10-29"
   
   日期格式設定可輸入help(“strtime”)查閱

檔案處理

• 自外部讀入一個文字檔非常容易, 只要一個指令
  

mydata=read.table("J:\\R\\Programs\\testdata.dat",header=TRUE, sep=",")

上例第一個參數為文字檔名, 第二個參數指明文字檔第一列為欄位名稱, 第三個參數指出欄位間分隔符號為逗點(,),其他可設定的參數非常多,可輸入指令help("read.table")查閱。

• 自R輸出一個文字檔只要改read為write即可
  

write.table(mydata,"J:\\R\\Programs\\testdata2.dat", sep=",",row.names = FALSE)

其中row.names = FALSE表示不輸出列號, 否則, R預設會輸出列號。

• 大檔案的處理

大檔一次讀入R, 會造成記憶體不足, 這時必須採取串流的方式, 一行行或一段段的讀, 才不會出問題, 指令如下

> con=file("j:/r/programs/testdata.dat","r") #開啟檔案

> readLines(con,1) #1表讀一列, 可一次讀取100列

> close(con) #關閉檔案