时间序列的简单R函数和例子.doc

一：时间序列的预处理1. scan从键盘或外部文件读入向量或列表。常用的形式为data-scan(file=“filepath&name”，skip=1,sep=“ ) 其它函数read.table， readLines ，write等。2. ts用来构造一个R中的时间序列对象。常用形式z - ts(data, start=c(1961, 1), frequency=12)jj- ts(jj, start=1960, frequency=4)3. plotR中的绘图工具。其外还有line,points等函数，屏幕分割函数par(mfrow=c(2,1)plot(x,y, type = p, col = red, lwd=10, main= 散点图 )4. acf计算样本自协方差函数，自相关函数以及偏相关系数。acf(x, lag.max = NULL, type = c(correlation, covariance, partial), plot = TRUE, na.action = na.fail, demean = TRUE, .)5. filter生成时间序列的线性滤波filter(x, filter, method = c(convolution, recursive), sides = 2)recursive 模式是自回归的，第0步系数默认为1， yi = xi + f1*yi-1 + . + fp*yi-pThe convolution filter is yi = f1*xi+o + . + fp*xi+o-(p-1)其中o 是偏移量，依赖sides的值。 6. lm线性回归函数lm(formula)没有截距的回归可用y x - 1 或 y 0 + x表示7. decompose利用滑动平均获取趋势项的分解方法，更复杂的有stl函数。decompose(x, type = c(additive, multiplicative), filter = NULL)8. Box.test 检验序列是否独立的函数Box.test(x, lag = 1, type = c(Box-Pierce, Ljung-Box)实例Example 1.1 jj = scan(jj.dat) # yes forward slash jj=ts(jj,start=1960, frequency=4) plot(jj, ylab=Quarterly Earnings per Share)Example 1.9 w = rnorm(500,0,1) # 500 N(0,1) variates v = filter(w, sides=2, rep(1,3)/3) # moving average par(mfrow=c(2,1) plot.ts(w) plot.ts(v)Example 1.10 w = rnorm(550,0,1) # 50 extra to avoid startup problems x = filter(w, filter=c(1,-.9), method=recursive) plot.ts(x51:550)Example 1.25 soi=scan(soi.dat) rec=scan(recruit.dat) par(mfrow=c(3,1) acf(soi, 50) acf(rec, 50) ccf(soi, rec, 50)Example 2.1 gtemp = scan(globtemp.dat) x = gtemp45:142 # use years 1990 to 1997 t = 1900:1997 fit=lm(xt) # regress x on t summary(fit) # regression output plot(t,x, type=o, xlab=year, ylab=temp deviation) abline(fit) # add regression line to the plotExample 2.2plot(decompose(AirPassengers)二：ARIMA建模相应函数9. arima.sim从arima模型得到模拟数据。arima.sim(n = 63, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488), sd = sqrt(0.1796) )10. ARMAacf， ARMAtoMA，acf2ARARMAacf利用参数计算理论自相关函数或偏自相关函数；ARMAtoMA利用模型参数估计wold系数；acf2AR利用自相关系数计算模型自回归模型系数。ARMAacf(ar = numeric(0), ma = numeric(0), lag.max = r, pacf = FALSE)ARMAtoMA(ar = numeric(0), ma = numeric(0), lag.max)acf2AR(Acf)11. ar ，ar.ols，ar.yw，ar.mle利用最小二乘,yule-walker,极大似然估计方法估计自回归模型的参数，默认利用AIC进行模型复杂度的选择。如果指定阶数，需要给定order.max的值。ar(x, aic = TRUE, order.max = NULL, method=c(yule-walker, burg, ols, mle, yw), na.action, series, .)ar.ols(x, aic = TRUE, order.max = NULL, na.action = na.fail, demean = TRUE, intercept = demean, series, .)12. arima利用数据拟合arima模型arima(USAccDeaths, order = c(0,1,1), seasonal = list(order=c(0,1,1)13. predict14. tsdiag给出时间序列模型的残差的诊断。fit - arima(lh, c(1,0,0)tsdiag(fit)实例Example 3.9 set.seed(5) ar2 = arima.sim(list(order = c(2,0,0), ar =c(1.5,-.75), n = 144) par(mfrow=c(3,1)plot.ts(ar2, axes=F); box(); axis(2) axis(1, seq(0,144,24) abline(v=seq(0,144,12), lty=dotted) acf = ARMAacf(ar=c(1.5,-.75), ma=0, 50) plot(acf, type=h, xlab=lag) abline(h=0)pacf = ARMAacf(ar=c(1.5,-.75), ma=0, 50, pacf=T)plot(pacf, type=h, xlab=lag) abline(h=0)Example 3.10 ARMAtoMA(ar=.9, ma=.5, 50) # for a list plot(ARMAtoMA(ar=.9, ma=.5, 50) # for a graphExample 3.16 rec = scan(recruit.dat) par(mfrow=c(2,1) acf(rec, 48) pacf(rec, 48) fit=ar.ols(rec,aic=F,order.max=2,demean=F,intercept=T) fit # estimates fit$asy.se # standard errorsExample 3.26 rec.yw = ar.yw(rec, order=2) rec.yw$x.mean rec.yw$ar sqrt(diag(rec.yw$asy.var.coef) rec.yw$var.pred rec.pr = predict(rec.yw, n.ahead=24) U = rec.pr$pred + rec.pr$se L = rec.pr$pred - rec.pr$se month = 360:453 plot(month, recmonth, type=o, xlim=c(360,480), ylab=recruits) lines(rec.pr$pred, col=red, type=o) lines(U, col=blue, lty=dashed) lines(L, col=blue, lty=dashed)Example 3.35 gnp96 = read.table(gnp96.dat) gnp = ts(gnp96,2, start=1947, frequency=4) par(mfrow=c(3,1)plot(gnp) acf(gnp, 50) gnpgr = diff(log(gnp) # growth rate plot.ts(gnpgr) par(mfrow=c(2,1) acf(gnpgr, 24) pacf(gnpgr, 24)# ARIMA fits: gnpgr.ar = arima(gnpgr, order = c(1, 0, 0) gnpgr.ma = arima(gnpgr, order = c(0, 0, 2)# to view the results: gnpgr.ar # potential problem here (see R Issue 1) gnpgr.ma ARMAtoMA(ar=.35, ma=0, 10) # prints psi-weightsExample 3.36 tsdiag(gnpgr.ma, gof.lag=20)Example 3.39 # AIC gnpgr.ma$aic # MA(2) gnpgr.ar$aic # AR(1)Example 3.43 prod=scan(prod.dat) par(mfrow=c(2,1) # (P)ACF of data acf(prod, 48) pacf(prod, 48) par(mfrow=c(2,1) # (P)ACF of d1 data acf(diff(prod), 48) pacf(diff(prod), 48) par(mfrow=c(2,1) # (P)ACF of d1-d12 data acf(diff(diff(prod),12), 48) pacf(diff(diff(prod),12), 48)# fit model (iii) prod.fit3 = arima(prod, order=c(1,1,1), seasonal=list(order=c(2,1,1), period=12) prod.fit3 # to view the results tsdiag(prod.fit3, gof.lag=48) # diagnostics# forecasts for the final model prod.pr = predict(prod.fit3, n.ahead=12) U = prod.pr$pred + 2*prod.pr$se L = prod.pr$pred - 2*prod.pr$se month=337:372 plot(month, prodmonth, type=o, xlim=c(337,384), ylim=c(100,180), ylab=Production) lines(prod.pr$pred, col=red, type=o) lines(U, col=blue, lty=dashed) lines(L, col=blue, lty=dashed) abline(v=372.5,lty=dotted)# Note: This can be done using acf2.R, sarima.R # and sarima.for.R (over here) as follows # (remember to source the code first): prod=scan(prod.dat) acf2(prod,48) acf2(diff(prod), 48) acf2(diff(diff(prod),12), 48)# fit model (iii) sarima(prod,1,1,1,2,1,1,12)# forecasts for the final model sarima.for(prod,12,1,1,1,2,1,1,12)#Example 4.9 soi = scan(soi.dat) rec = scan(recruit.dat) par(mfrow=c(2,1) soi.per = spec.pgram(soi, taper=0, log=no) abline(v=1/12, lty=dotted) abline(v=1/48, lty=dotted) rec.per = spec.pgram(rec, taper=0, log=no) abline(v=1/12, lty=dotted) abline(v=1/48, lty=dotted) soi.per$spec40 # soi pgram at freq 1/12 = 40/480 soi.per$spec10 # soi pgr