機關根檢驗urdf_R語言時間序列函數整理[轉]]

文檔1：《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")

[email protected]  #統計量

[email protected]      #臨界值

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")

[email protected]  #統計量

[email protected]      #臨界值

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

}

#...