R語言時間序列函數整理__函數

【包】
library(zoo)            #時間格式預先處理
library(xts)            #同上
library(timeSeires)      #同上
library(urca)           #進行單位根檢驗
library(tseries)         #arma模型
library(fUnitRoots)     #進行單位根檢驗
library(FinTS)         #調用其中的自迴歸檢驗函數
library(fGarch)        #GARCH模型
library(nlme)          #調用其中的gls函數
library(fArma)        #進行擬合和檢驗

【基本函數】
數學函數
abs，sqrt：絕對值，平方根 log, log10, log2 , exp：對數與指數函數 sin，cos，tan，asin，acos，atan，atan2：三角函數 sinh，cosh，tanh，asinh，acosh，atanh：雙曲函數 
簡單統計量
sum, mean, var, sd, min, max, range, median, IQR（四分位間距）等為統計量，sort，order，rank與排序有關，其它還有ave，fivenum，mad，quantile，stem等。


【資料處理】
#具體說明見文檔1
#轉成時間序列類型
x = rnorm(2)
charvec = c(“2010-01-01”,”2010-02-01”)
zoo(x,as.Date(charvec))     #包zoo
xts(x, as.Date(charvec))     #包xts
timeSeries(x,as.Date(charvec))  #包timeSeries
#規則的時間序列，資料在規定的時間間隔內出現
tm = ts(x,start = c(2010,1), frequency=12 )  #12為按月份，4為按季度，1為按年度
zm = zooreg(x,start = c(2010,1), frequency=12 )  #包zoo
xm = as.xts(tm)     #包xts
sm = as.timeSeries(tm) #包timeSeries
#判斷是否為規則時間序列
is.regular(x)

#排序
zoo()和xts()會強制變換為正序（按照時間名稱）
timeSeries不會強制排序；其結果可以根據sort函數排序，也可以採用rev()函數進行逆序；參數recordIDs，可以給每個元素（行）標記一個ID，從而可以找回原來的順序

#預設的時間有重複的時間點時
zoo會報錯
xts按照升序排列
timeSeries把重複部分放置在尾部； 

#行合并和列合并
#都是按照列名進行合并，列名不同的部分用NA代替
cbind()
rbind()
merge() 列合并

#取子集
xts()預設將向量做成了矩陣；其他與常規向量或者矩陣沒有差別

#缺失值處理
na.omit(x) 
x[is.na(x)] = 0
x[is.na(x)] = mean(x,na.rm=TRUE)
x[is.na(x)] = median(x,na.rm=TRUE)
na.approx(x)  #對缺失值進行線性插值
na.spline(x)   #對缺失值進行樣條插值
na.locf(x)     #末次觀測值結轉法
na.trim(x, sides=”left” )  #去掉最後一個缺失值
#對timeSreies資料
na.omit(x, “ir” )  #去掉首末位置的缺失值
na.omit(x, “iz” )  #用替換首末位置的缺失值
na.omit(x, “ie” )  #對首末位置的缺失值進行插值
na.omit(x, method=“ie”, interp= c(“before”,”linear”,”after”) ) #可以選擇插值方法，before末次觀測值法，after下次觀測結轉法

as.contiguous(x)  #返回x中最長的連續無缺失值的序列片段，如果有兩個等長的序列片段，則返回第一個。

#時間序列資料的顯示
#zoo和xts都只能按照原來的格式顯示，timeSeries可以設定顯示格式
print(x, format= “%m/%d/%y %H:%M”)  #%m表示月，%d表示天，%y表示年，%H表示時，%M表示分鐘，%A表示星期，%j表示天的序號
#timeSeries也可以按照ts的格式顯示
print(x, style=”ts”)
print(x, style=”ts”, by=”quarter”)

【圖形展示】
plot.zoo(x)
plot.xts(x)
plot.zoo(x, plot.type=”single”) #支援多個時間序列資料在一個圖中展示
plot(x, plot.type=”single”) #支援多個時間序列資料在一個圖中展示，僅對xts不行

【基本統計運算】
1、自相關係數、偏自相關係數等
例題2.1 
d=scan("sha.csv")
sha=ts(d,start=1964,freq=1)
plot.ts(sha)   #繪製時序圖
acf(sha,22)   #繪製自相關圖，滯後期數22
pacf(sha,22)  #繪製偏自相關圖，滯後期數22
corr=acf(sha,22)   #儲存相關係數
cov=acf(sha,22,type = "covariance")   #儲存共變數

2、同時繪製兩組資料的時序圖
d=read.csv("double.csv",header=F)
double=ts(d,start=1964,freq=1)
plot(double, plot.type = "multiple")   #兩組資料兩個圖
plot(double, plot.type = "single")     #兩組資料一個圖 
plot(double, plot.type = "single",col=c("red","green"),lty=c(1,2)) #設定每組資料圖的顏色、曲線類型)

3、純隨機性檢驗
例題2.3續
d=scan("temp.csv")
temp=ts(d,freq=1,start=c(1949))
Box.test(temp, type="Ljung-Box",lag=6)

4、差分運算和滯後運算
diff
lag

5、類比ARIMA模型的結果
arima.sim(n = 100, list(ar = 0.8))
plot.ts(arima.sim(n = 100, list(ar = 0.8)))   #會隨機產生一個包含100個隨機數的時序圖
plot.ts(arima.sim(n = 100, list(ar = -1.1)))   #非平穩，無法得到時序圖。
plot.ts(arima.sim(n = 100, list(ar = c(1,-0.5))))
plot.ts(arima.sim(n = 100, list(ar = c(1,0.5))))
arima.sim(n = 1000, list(ar = 0.5, ma = -0.8))
acf(arima.sim(n = 1000, list(ar = 0.5, ma = -0.8)),20)
pacf(arima.sim(n = 1000, list(ar = 0.5, ma = -0.8)),20)

【單位根檢驗】
#方法1
b=ts(read.csv("6_1.csv",header=T)) 
x=b[,1]
y=b[,1]
summary(ur.df(x,type="trend",selectlags="AIC"))
#方法2：單位根檢驗更好的函數，加了畫圖的功能 
library(fUnitRoots)
urdfTest(x)
#方法3：ADF檢驗的一個自編函數
library(urca)
#...
ur.df.01=function(x,lags=8){    
  #將三種ADF檢驗形式匯總的函數（結果和EVIEWS不一致）
  res=matrix(0,5,3)
  colnames(res)=c("無","含常數項","含常數項和趨勢項")
  rownames(res)=c("tau統計量","1%臨界值","5%臨界值",
                  "10%臨界值","是否穩定（1/0）")
  types=c("none","drift","trend")
  for(i in 1:3){
    x.adf=ur.df(x,type=types[i],lags=lags,selectlags="AIC")
    x.adf.1=x.adf@teststat  #統計量
    x.adf.2=x.adf@cval      #臨界值
    res[1,i]  =x.adf.1[1]
    res[2:4,i]=x.adf.2[1,]
    res[5,i]=if( abs(res[1,i]) > abs(res[3,i]) ) 1 else 0
  }
  return(res)
}
#...
ur.df.01(x)              #對原序列進行判斷

【一般的ARIMA模型】
d=scan("a1.5.txt")               #匯入資料
prop=ts(d,start=1950,freq=1)      #轉化為時間序列資料
plot(prop)                   #作時序圖
acf(prop,12)                 #作自相關圖，拖尾
pacf(prop,12)                #作偏自相關圖，1階截尾
Box.test(prop, type="Ljung-Box",lag=6)  
#純隨機性檢驗,p值小於5%，序列為非白色雜訊
Box.test(prop, type="Ljung-Box",lag=12)
( m1=arima(prop, order = c(1,0,0),method="ML") )    #用AR(1)模型擬合，如參數method="CSS"，估計方法為條件最小二乘法，用條件最小二乘法時，不顯示AIC。
( m2=arima(prop, order = c(1,0,0),method="ML", include.mean = F) ) #用AR(1)模型擬合，不含截距項。
tsdiag(m1)  #對估計進行診斷，判斷殘差是否為白色雜訊
summary(m1)
r=m1$residuals  #用r來儲存殘差
Box.test(r,type="Ljung-Box",lag=6, fitdf=1)#對殘差進行純隨機性檢驗，fitdf表示殘差減少的自由度
AutocorTest(m1$resid)                    #載入FinTS包，進行自相關檢驗
prop.fore = predict(m1, n.ahead =5)  #將未來5期預測值儲存在prop.fore變數中
U = prop.fore$pred + 1.96* prop.fore$se  #會自動產生方差
L = prop.fore$pred – 1.96* prop.fore$se   #算出95%信賴區間
ts.plot(prop, prop.fore$pred, col=1:2)      #作時序圖，含預測。
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")#在時序圖中作出95%信賴區間

——說明：運行命令arima(prop, order = c(1,0,0),method="ML")之後，顯示：
Call:
arima(x = prop, order = c(1, 0, 0), method = "ML")
Coefficients:
         ar1    intercept
      0.6914    81.5509
s.e.   0.0989     1.7453
sigma^2 estimated as 15.51:  log likelihood = -137.02,  aic = 280.05
注 意：intercept下面的81.5509是均值，而不是截距。雖然intercept是截距的意思，這裡如果用mean會更好。（the mean and the intercept are the same only when there is no AR term，均值和截距是相同的，只有在沒有AR項的時候）
如果想得到截距，利用公式計算。int=(1-0.6914)*81.5509= 25.16661。

——說明：Box.test(r,type="Ljung-Box",lag=6，fitdf=1)
fitdf表示p+q，number of degrees of freedom to be subtracted if x is a series of residuals，當檢驗的序列是殘差到時候，需要加上命令fitdf，表示減去的自由度。
運行Box.test(r,type="Ljung-Box",lag=6,fitdf=1)後，顯示的結果：
Box.test(r,type="Ljung-Box",lag=6,fitdf=1)
        Box-Ljung test
data:  r 
X-squared = 5.8661, df = 5, p-value = 0.3195
“df = 5”表示自由度為5，由於參數lag=6，所以是滯後6期的檢驗。

#另一個參數估計與檢驗的方法（載入fArma程式包）
ue=ts(scan("unemployment.txt"),start=1962,f=4) #讀取資料
due=diff(ue)
ddue=diff(due,lag=4)
fit2=armaFit(~arima(4,0,0),include.mean=F,data=ddue,method="ML")  #另一種擬合函數
summary(fit2)
fit3=armaFit(~arima(4,0,0),data=ddue,transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="CSS")
summary(fit3)

【一些特殊的模型】
#固定某些係數的值
arima(dw,order=c(4,0,0),fixed=c(NA,0,0,NA,0),method="CSS")

#乘積季節模型
wue=ts(scan("wue.txt"),start=1948,f=12)
arima(wue,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12),include.mean=F,method