r時間序列分析 為什麼定階數,如何定,如何判斷
r時間序列分析工具
xts包 xts(x=null,order.by=index(x),…) coredata() xts資料子集
ohlc資料格式
quantmod包
ttr包
自回歸模型(ar) 跟以前時刻有關和當前隨機游動有關
ar(p)的性質
平穩性要求：
ar(2)特徵根的模都小於1
ar(p)的定階 赤池資訊準則（aic） 施瓦茨—貝葉斯資訊準則（bic）
aic bic最小值的地方定階
白雜訊時是嚴平穩序列
平穩時間序列
# ar models 自回歸模型ar(p) p階
da=read.table("q-gnp4710.txt",header=t)
head(da)
g=da$value
lg=log(g) #取對數
gnp=diff(lg) #lg的差分,返回適當的滯後，迭代的差異
dim(da) #重新設定或設定物件的維數
tdx=c(1:253)/4+1947 # create the time index
par(mfcol=c(2,1))
plot(tdx,g,xlab='year',ylab='gnp',type='l')
plot(tdx[2:253],gnp,type='l',xlab='year',ylab='growth')
acf(gnp,lag=12) #計算互相關或兩個單變數序列的協方差
pacf(gnp,lag=12) # compute pacf
m1=arima(gnp,order=c(3,0,0)) #擬合arima模型的單變數時間序列 order=c(3,0,0)第乙個是ar階數第二個是差分階數第三個是ma模型階數
m1tsdiag(m1,gof=12) # model checking discussed later繪製時間序列
p1=c(1,-m1$coef[1:3]) # set-up the polynomial設定多項式
r1=polyroot(p1) # solve the polynomial equation求解多項式方程
r1mod(r1) #絕對值 (模)
k=2*pi/acos(1.616116/1.832674) # compute length of the period 計算週期長度
kbox.test(m1$residuals,lag=4,type='ljung') #模型檢驗
1-pchisq(0.74926,1)
forecast(m1,1) #**第一步
forecast(m1,2) #**第二步
滑動平均模型(ma) 與上幾個隨機游動有關
da=read.table("m-ibm3dx2608.txt",header=t)
head(da)
ew=da$ewrtn
plot(da$data,ew,type="1")
#定階 取不為0的階不要太高(複雜)不要太低(不准)非0說明相關
acf(ew,lag=20)
#模型擬合
m1=arima(ew,order=c(0,0,9))
m1m1=arima(ew,order=c(0,0,9),fixed=c(na,0,na,0,0,0,0,0,na,na))#只擬合相關係數不為0的項 aic越小越好
m1sqrt(0.005097) #計算隨機游動標準差
#模型檢驗
tsdiag(m1,gof=12) #檢測殘差是否為白雜訊，殘差為白雜訊好
box.test(m1$residuals,lag=12,type='ljung') # model checking
pv=1-pchisq(17.6,9) # compute p-value after adjusting the d.f.
pv#**
m1=arima(ew[1:986],order=c(0,0,9),fixed=c(na,0,na,0,0,0,0,0,na,na))#只擬合相關係數不為0的項 aic越小越好
m1predict(m1,10)
pre=c(predict(m1,10)$pred)
data.frame(ew[987:996],pre,ew[987:996]-pre) #**值與真實值比較
自回歸滑動平均模型(arma)
ar(p)模型與ma(q)模型的結合——arma(p,q)
arma(p,q)定階
推廣的自相關函式（eacf） 非零
#構建arma模型
da=read.table("m-ibm3dx2608.txt",header=t)
head(da)
ew=da$ewrtn
acf(ew)
pacf(ew)
#根據資訊準則定階
library("forecast")
auto.arima(ew)
#根據eacf定階
eacf(ew)
#模型擬合
m2=arima(ew,order=c(2,0,2))
m2m3=arima(ew,order=c(0,0,1))
m3#模型檢驗
tsdiag(m2,gof=12) #檢測殘差是否為白雜訊
box.test(m2$residuals,lag=12,type='ljung') # model checking
tsdiag(m3,gof=12) #檢測殘差是否為白雜訊
box.test(m3$residuals,lag=12,type='ljung') # model checking
#**predict(m2,1)
predict(m3,1)
平穩性檢驗
adf單位根檢驗library("funitroots")
adftest()
時間序列 確定性趨勢成分+平穩的零均值誤差成分
確定性因素分解
長期趨勢變動
迴圈變動
季節性變動
隨機波動
確定性時間序列分析
– 平滑**法（平穩時間序列）
簡單移動平均法sma
k期移動平均法
k期中心移動平均
– 指數平滑法
– 趨勢**法
線性擬合
曲線擬合
– 分解分析法
加權移動平滑（wma）
– k期移動平均
– 指數加權移動平滑（ewma）
簡單移動平均法sma
sma.cal <- function( ts )
加權移動平均法(時間越近,權重越高)
wma.cal <- function( ts, weight )
k期移動平均法(簡單移動平均法的優化,取最近k期的均值)
kwma.cal <- function( ts, k=20 )
指數平滑法 平滑係數α的選擇(時間越遠權重小)
ewma.cal <- function( ts, a=0.8 )
aa=c(1:10)
n=length(aa)
a=0.8
pred=0
pred[1:3]<-aa[1:3]
for(i in 4:length(aa))#pred為**值，其中前三個資料為真實資料
p=a*aa[n]+(1-a)*pred[n]
季節均值法
余弦趨勢法
季節指數水平法
季節指數趨勢法
全年比率平均法
步驟： 1. 求各年平均數
2. 求各月對該年平均數比率
3. 計算季節比率
4. **
非平穩時間序列
平穩性檢驗 adf單位根檢驗 adftest()
library(funitroots)
adftest(oil.price,lags=12,type="c") #得出p值 取到這個樣本的概率是多少 p值即概率，反映某一事件發生的可能性大小
引數lags(滯後項)的設定——根據ar模型的階數估算而確定
data("oil.price") #tsa包裡有
ar(oil.price,method = 'mle')$order #得出階數
識別趨勢非平穩和差分非平穩
趨勢非平穩
差分非平穩
library(funitroots)
da=read.table("d-sp55008.txt",header=t)
head(da)
sp5=log(da[,7])
m2=ar(diff(sp5),method='mle') # based on aic
m2$order
adftest(sp5,lags=2,type=("ct"))
adftest(sp5,lags=15,type=("ct")) # based on pacf
非平穩ar(1)模型 係數<1
data(explode.s)
plot(explode.s,ylab=expression(y[t]),type="o")
acf(explode.s)
pacf(explode.s)
對於隨機游動 不易過度差分
一次差分後 是乙個白雜訊過程 在這裡就可以停止了不易過度差分
兩次差分後 ma(1)模型 過差分效果不好
arima(p,d,q)模型
如果乙個時間序列{
R語言 時間序列分析
時間序列分析 1.對時間序列的描述 2.利用前面的結果進行 ts是時間序列的英文簡稱 可以使用sys.date 函式檢視當前系統的時間 用seq函式創造連續的時間點 ts函式生成時間序列，可以很方便的將向量轉化成時間序列 egsales round runif 50,min 50,max 100 s...
《應用時間序列分析 R軟體陪同》 導讀
前言 首先,一些教材偏重於數學理論和推導.作者多為數學出身,他們習慣於數學的嚴格性和匯出精確而又漂亮的數學結論.這些書適用於那些願意為時間序列的數學理論研究做出貢獻的讀者.其次,國內教材中一元時間序列往往佔絕大部分篇幅,而且包含在各種數學假定下的各種定理和結果.這是因為一元時間序列的數學描述確實很漂...
R 簡述時間序列（ARIMA）
時間序列是指將同一統計指標的數值按時間先後的順序排列而成的數列，主要目的是用歷史資料 未來。時間序列一般分為連續性和離散型，下面只講離散型。我這裡主要講的是r中arima模型 單積自回歸移動平均 其中ar是自回歸 ma是移動平均，具體推導很複雜就不詳述了。白雜訊 均值，方差不隨時間變化的一種平穩的隨...