【包】
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="CSS")
#擬合自回歸模型，因變量關於時間的回歸模型
eg1=ts(scan("582.txt"))
ts.plot(eg1)
fit.gls=gls(eg1~-1+time(eg1), correlation=corARMA(p=1), method="ML") #看nlme包
summary(fit.gls2)
#或
fit=arima(eg1,c(1,0,0),xreg=time(eg1),include.mean=F,method="ML")
AutocorTest(fit$resid) #殘差白噪聲檢驗
#延遲因變量回歸模型
leg1=lag(eg1,-1)
y=cbind(eg1,leg1)
fit=arima(y[,1],c(0,0,0),xreg=y[,2],include.mean=F)
#擬合GARCH模型
library(tseries)
library(fGarch)
library(FinTS)
a=ts(scan("583.txt"))
ts.plot(a)
fit=lm(a~-1+time(a))
r=resid(fit)
summary(fit)
pacf(r^2)
acf(r)
acf(r^2)
AutocorTest(r) #殘差是否存在序列相關
ArchTest(r) #是否存在ARCH效應
fit1=garchFit(~arma(2,0)+garch(1,1), data=r, algorithm="nlminb+nm",
trace=F, include.mean=F)
summary(fit1)
#協整檢驗
fit=arima(b[,2],xreg=b[,1],method="CSS")
r=resid(fit)
summary(ur.df(r,type="drift",lag=1))
Box.test(r,lag=6,fitdf=1)
【自動運行的自編函數】
acf.3(x) #同時繪製3個相關圖，acf函數的擴展
ur.df.01(x) #進行單位根檢驗，得到更加舒服的結果
tsdiag2(x) #返回x的
arma.choose(x,ari=3,mai=3) #選擇合適的AR和MA，基於包tseries的arma函數
#########################附屬自編函數
#...
acf.3=function(x,lag.max=10,…){
ol=par(mfrow=c(3,1),mar=c(2,4,1,1))
acf(x,lag.max=lag.max,type="correlation")
acf(x,lag.max= lag.max,type="covariance")
acf(x,lag.max= lag.max,type="partial")
par(ol)
}
#...
#...類似於tsgiag函數的擴展
tsdiag2=function(xx.model,fitdf=0,testlag=10){
t1=xx.arma$residuals
t2=acf(na.omit(t1),plot=F)
t3=sapply(1:testlag,
function(x,r,fitdf){
Box.test(r,type="Ljung-Box",lag=x, fitdf=fitdf)
},
r=t1,fitdf=fitdf)
par(mfrow=c(3,1))
plot(t1,type="b",ylab="",main="殘差走勢")
lines(c(0,length(t1)*2),c(0,0),col=2,lty=2)
plot(t2,type="h",ylab="ACF",main="殘差的自相關係數")
plot(do.call("c",t3[3,]),type="p",ylab="P-value",pch=16,col=4,
ylim=c(0,1),main="殘差的Ljung-Box檢驗")
lines(c(0,attr(t1,"tsp")[2]),c(0.05,0.05),lty=2,col=2)
}
#...
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( !is.nan(res[1,i]) & abs(res[1,i]) > abs(res[3,i]) ) 1 else 0
}
return(res)
}
#...
#...
arma.choose.02=function(x){
#二進位進位運算，以矩陣形式,x=c(0,1,0,1,...)
n=length(x)
if( all(!as.logical(x-rep(1,n))) ) stop("已不能再加1!")
x[1]=x[1]+1
for(i in 1:(n-1)) if(x[i]>1){ x[i]=0;x[i+1]=x[i+1]+1 }
return(x)
}
arma.choose.01=function(ti){
#把ti變換成所有可能的ti個0或1的組合
if(ti<0) stop("ti要大於0！")
if(ti==0) return(0)
if(ti%%1!=0) stop("ti要整數！")
res=matrix(0,2^ti,ti)
for(i in 2:2^ti) res[i,]=arma.choose.02(res[i-1,])
return(res)
}
arma.choose.03=function(t0){
gsub(", ",".",toString(t0,sep=""))
}
arma.choose.04=function(i,ari,tti){
#ari是最大滯後期,tti由ari生成
ar.lag=((1:ari)*tti[i,])
ar.lag=ar.lag[ar.lag!=0]
ar.lag
}
arma.choose=function(x,ari=3,mai=3,...){
tti=arma.choose.01(ari)
ttj=arma.choose.01(mai)
ti=2^ari;tj=2^mai
res.aic=matrix(Inf,ti,tj) #保存所有組合的AIC
rownames(res.aic)=paste("AR",apply(tti,1,arma.choose.03),sep=".")
colnames(res.aic)=paste("MA",apply(ttj,1,arma.choose.03),sep=".")
res.rss=matrix(Inf,ti,tj) #保存所有組合的RSS
rownames(res.rss)=paste("AR",apply(tti,1,arma.choose.03),sep=".")
colnames(res.rss)=paste("MA",apply(ttj,1,arma.choose.03),sep=".")
for(i in 2:ti){
j=1
ar.lag=arma.choose.04(i,ari,tti)
x.arma=arma(x,lag=list(ar=ar.lag),...)
ss=summary(x.arma)
res.aic[i,j]=ss$aic
res.rss[i,j]=sum(ss$residuals^2)
}
for(j in 2:tj){
i=1
ma.lag=arma.choose.04(j,mai,ttj)
x.arma=arma(x,lag=list(ma=ma.lag),...)
ss=summary(x.arma)
res.aic[i,j]=ss$aic
res.rss[i,j]=sum(ss$residuals^2)
}
for(i in 2:ti){for(j in 2:tj){
ar.lag=arma.choose.04(i,ari,tti)
ma.lag=arma.choose.04(j,mai,ttj)
x.arma=arma(x,lag=list(ar=ar.lag,ma=ma.lag),...)
ss=summary(x.arma)
res.aic[i,j]=ss$aic
res.rss[i,j]=sum(ss$residuals^2)
}}
res=list()
res[["tt.ar"]]=tti
res[["tt.ma"]]=ttj
temp1=which.min(res.aic) #找到最小的位置，把res.aic當做按列排的向量
temp2=temp1 %% ti #ti是行數，取餘以後就是（temp2）行號
#AR可以直接被arma調用，MA同理
res[["AR"]]=if(temp2==0) arma.choose.04(ti,ari,tti) else arma.choose.04(temp2,ari,tti)
res[["MA"]]=arma.choose.04( ceiling( temp1 / ti ), mai,ttj)
res[["aic"]]=res.aic
res[["rss"]]=res.rss
res
}
#