2008年11月8日 星期六
時間數列 (二) – 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,變異數為σ2,且εi及εj無關聯
之後我們用誤差的最小平方法(Least Square),求出a及b,進而得到
- 總平方和(SST) = 殘差平方和(SSE) + 迴歸平方和(SSR)
- R2 = SSR / SST,R2是檢定x與y的相關係數
- 檢定一元迴歸可用自由度=n-2的t分配,多元迴歸則用F分配
F = (SSR / (k -1)) / (SSE / (n – k))
其中k-1為自變數個數
時間數列與迴歸分析不同點就是假設誤差不是完全隨機的, 而是自我相關的,即前期的誤差會影響後期的誤差
εt = ρεt -1 + at
其中ρ絕對值<1,at稱為干擾項(Disturbance Term),服從常態分數,平均數為0,變異數為σ2
- 檢定ρ是否為0,一階(誤差只與前一期有關)可用Dubin-Watson統計量檢定,多階可用ARIMA統計量。
- ARIMA: Auto Regression Integrated Moving Average,為Box/Jenkins提出
此外,若考慮不是σ2固定常數,也就是異質變異數(Heteroscendasticity),可用加權最小平方法處理,這就是ARCH或GARCH模型。
註:此文主要是節錄自林茂文先生所注之時間數列之分析與預測一書。
2008年11月1日 星期六
迴歸(Regression)
就從迴歸作起吧.
- 載入測試資料, 也可以自己輸入
> load("j:/r/programs/testdata/usingr.rdata") - 使用 moths 變數, 先輸入moths, 查看資料
- 以moths$meters 為應變數, moths$A、moths$P 為自變數
> lm(moths$meters ~ moths$A+moths$P) - 得到簡單迴歸的截距(Interception)及斜率(Slope)
Coefficients:(Intercept) moths$A moths$P
45.1017 0.3712 -0.6005
R 變數類型
- 向量: 為預設類型, 即一維陣列
- Data frame: 具有欄位名稱的二維陣列
- List: 可儲存不同類型的物件集合, 例如summary(), 要轉成向量, 可用unlist()
- 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) - S4: 類似S3, 有點複雜, 以後再說明
- 日期變數: 要使用日期函數處理, 必須先把變數轉為日期
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預設會輸出列號。
- 大檔案的處理
> con=file("j:/r/programs/testdata.dat","r") #開啟檔案
> readLines(con,1) #1表讀一列, 可一次讀取100列
> close(con) #關閉檔案
R 常用指令
- 變數指定
> a=1 - 變數預設資料型態為向量, 亦稱為集合(collection), 故使用c
> x = c(10.4, 5.6, 3.1, 6.4, 21.7) - 特殊變數指定方式
> x = 1:12 # x=1,2,3,4, ..., 12
> a = seq(-5,5,by=.2) # -5~5間隔0.2產生一個值 - 計算結果
> a * 2 - 列出所有變數
> ls() - 列出所有變數及資料型態(data type)
> ls.str() - 註解(Comment): # 開頭
- 幫變數加註解
> comment(x1) = c("This data is from experiment #0234", "Jun 5, 1998") - 移除變數 c, d
> rm(c,d) - 移除所有變數
> rm(list=ls()) - 查看指令說明
> help(“rm”) - 搜尋相關說明
> help.search("Time Series") - 查看套件說明
> library(help="stats") - 二維資料型態 data frame, 它具有欄位名稱, 如下例, 欄位名稱分別為f1、f2, 資料有三列
> a = data.frame(f1=c(1,2,3),f2=c(2,4,6)) - 取其中一欄的值
> a$f1 - 兩個向量合併為一data frame
> x=cbind(a,b) # 若 a, b 元素個數不等, R會自動重複元素, 已補齊至最大元素個數
> x=rbind(a,b) # a, b 元素已列方式合併 - 以特定欄位比對合併兩data frame為一data frame
#以兩data frame 的日期比對, 合併
> x2=merge(mydata,mydata2, by.x="Date", by.y="Date")
判斷資料型態是否為矩陣
> is.matrix(b) - 顯示資料型態
> class(b) - 矩陣轉置(transpose)
> t(b) - 矩陣相乘, 須用%*%, 若用*, 系統會將相同位置的兩數相乘, 而非矩陣相乘
> t(b)%*%b - 引入一段R程式, 假設你寫好一段公用函數, 你可以以下指令執行, 而不必每次重打, 或複製/貼上方式執行
source("c:/1.r")
R 的環境簡介
個人認為 R 的強大功能有兩項:
- 專為數學/統計設計的程式語言, 變數不必宣告, 直接使用, 函數庫(Library)不必指定, 只需載入即可, 向量、矩陣、微積分、統計分配函數及轉換樣樣不缺, 繪圖能力超強。
- 套件(Package)特多, 官方網站數量就超過百個, 也可以自製套件。
當然, 也不是沒有缺點, 套件說明稍嫌簡略, 不過這是Open Source的普遍現象吧.
在程式撰寫前, 應先了解開發環境:
- 不知道怎麼開始, 輸入help.start(), 即可連到各種使用手冊的首頁。
- 設定現行目錄: 執行選單 【檔案】>【變更現行目錄】或直接下指令如下, 設定程式存放的目錄, 目錄分隔符號(\)須改為斜線(/)或使用雙倒斜線(\\)。
setwd("J:\\R\\Programs") - 載入套件: 載入新套件前, 先下指令列出目前已載入的套件
library()
輸入下列指令載入套件
install.packages("graphics")
之後選擇下載的地點即可開始下載,並同時載入至目前環境中, 以上指令也可以執行選單, 達成相同效果。 - 儲存工作空間(Work Space): 輸入下列指令或執行選單 【檔案】>【儲存工作空間】, 將之前的設定儲存到一個組態檔(*.rdata)中。
save.image("J:\\R\\Programs\\My.rdata") - 載入工作空間(Work Space): 下次執行R時, 可先輸入下列指令或執行選單 【檔案】>【載入工作空間】, 載入上次的設定,
R的通則:
- 指令大小寫有區別
- 函數一律要加()
R 的初次邂逅
以下就簡單說明R的安裝與使用:
- 安裝: R的安裝很簡單, 可以在Windows或Linux安裝, 至 http://cran.r-project.org/ 下載並執行, 就搞定了。
- 啟動: 從Windows選單, 執行【開始】>【程式集】>【R】>【R 2.7.0】。
- 使用: 出現畫面如下, 在提示符號後輸入指令即可。例如
a=1
b=2
a+b - 離開: 輸入q(), 按ENTER即可。