R语言时间序列作业

1、2016 年第二学期时间序列分析及应用 R 语言课后作业 第三章 趋势3.4(a) data(hours);plot(hours,ylab='Monthly Hours',type='o')画出时间序列图(b) data(hours);plot(hours,ylab='Monthly Hours',type='l')type='o' 表示每个数据点都叠加在曲线上； type='b' 表示在曲线上叠加数据点，但是该数 据点附近是断开的； type='l' 表示只显示各数据点之间的连接线

2、段； type='p' 只想显示数据 点。points(y=hours,x=time(hours),pch=as.vector(season(hours).5.0M.5.04.0.004.5.0MJD1983 1984 198519861987Time3.10(a) data(hours);hours.lm=lm(hourstime(hours)+I(time(hours)2);summary(hours.lm)用最小二乘法拟合二次趋势，结果显示如下：Call: lm(formula = hours time(hours) + I(time(hours)2)Residuals:

3、Min 1Q Median 3Q Max-1.00603 -0.25431 -0.02267 0.22884 0.98358Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept)-5.122e+05 1.155e+05 -4.433 4.28e-05 *time(hours)5.159e+02 1.164e+02 4.431 4.31e -05 *I(time(hours)2)-1.299e -01 2.933e-02 -4.428 4.35e-05 *Signif. codes: 0 * ' 0.001 * &#

4、39; 0.01 * ' 0.05 .' 0.1 ' 1Residual standard error: 0.423 on 57 degrees of freedomMultiple R -squared: 0.5921, Adjusted R-squared: 0.5778F-statistic: 41.37 on 2 and 57 DF, p-value: 7.97e -12(b) plot(y=rstudent(hours.lm),x=as.vector(time(hours),type='l',ylab='Standardized Res

5、iduals') points(y=rstudent(hours.lm),x=as.vector(time(hours),pch=as.vector(season(hours)标准残差的时间序列，应用月度绘图标志。(为了更容易识别季节性)带季节性图标的的残差 -时间图slaudiseRdezidradnatS19831984 1985 19861987as.vector(time(hours)(c) runs(rstudent(hours.lm) 对标准差进行游程检验 $pvalue1 0.00012$observed.runs1 16$expected.runs1 30.96667

6、$n11 31 $n21 29$k1 0结果解释： P 值为 0.00012，表明非随机性是合理的。(d) acf(rstudent(hours.lm)标准残差的样本自相关函数季节均值模型残差的样本自相关系数Series rstudent(hours.lm)FCA10ga(e) qqnorm(rstudent(hours.lm);qqline(rstudent(hours.lm)(QQ 图)Normal Q-Q Plot-2 -1 0 1 2Theoretical Quantiles正态性可以通过正态得分或者分位数-分位数( QQ)图来检验。此处的直线型图形支持了该模型中随机项是正态分布的假设

7、。hist(rstudent(hours.lm),xlab='Standardized Residuals')标准残差的直方图(季节均值模型的标准残差直方图)Histogram of rstudent(hours.lm)Standardized Residualsshapiro.test(rstudent(hours.lm)正态性检验( Shapiro-Wilk 检验)本质是：计算残差与相应的正态分位数之间的相关系数。 相关性越小，就越有理由否定正态性。Shapiro-Wilk normality testdata: rstudent(hours.lm)W = 0.99385,

8、 p-value = 0.9909根据上面的检验结果，我们不能拒绝模型的随机项是正态分布的假设。第四章 平稳时间序列模型124.4第五章 非平稳时间序列模型5.1(a) ARMA (2,1) p=2,q=1,参数值和1=1 2=-0.25 1=0.1(b) IMA(2,0) p=2,d=1,q=0, 参数值和(c) ARMA(2,2) p=2,q=2, 参数值和1=0.52=-0.5 1=0.5 2=-0.255.7(a) A:AR(2) 1=0.92=0.09B:IMA(1,1) 1=0.1(b) 一个是固定的一个是不固定的。5.11(a) data(winnebago);win.graph

9、(width=6.5,height=3,pointsize=8)plot(winnebago,type='o',ylab='Winnebago Monthly Sales')时间序列图selaSylhtnoMogaben表明公司的休闲车的销量在逐渐增加。(b) plot(log(winnebago),type='o',ylab='Log(Monthly Sales)')取对数之后的时间序列图仍然呈现增加的趋势，但是比没有取对数之前增加的缓慢一些。(c) percentage=na.omit(winnebago -zlag(winn

10、ebago)/zlag(winnebago) win.graph(width=3,height=3,pointsize=8) plot(x=diff(log(winnebago) -1,y=percentage -1,ylab='PercentageChange',xlab='Difference of Logs')cor(diff(log(winnebago) -1,percentage -1)1 0.9646886Difference of Logs5.1 0.1 5.0 0.0 5.0egnahCegatnecreP结果显示： 0.96 认为一致。第六章

11、模型识别6.29(a) set.seed(762534);series=arima.sim(n=60,list(ar=0.4,ma=0.6)phi=0.4;theta=0.6;ACF=ARMAacf(ar=phi,ma= -theta,lag.max=10)plot(y=ACF -1,x=1:10,xlab='Lag',ylab='ACF',type='h',ylim=c( -.2,.2);abline(h=0)2.01.00.01.02.0FCA(b) acf(series)Series seriesLag第八章 模型诊断8.7(a) data

12、(hare);model=arima(sqrt(hare),order=c(3,0,0)win.graph(width=6.5,height=3,pointsize=8);acf(rstandard(model)残差的样本自相关图Series rstandard(model)1412101.01.03.0FCA(b) LB.test(model,lag=9)Box-Ljung testdata: residuals from modelX-squared = 6.2475, df = 6, p-value = 0.396 JB统计量结果表明不拒绝误差项的独立性。(c) 对残差进行检验。runs

13、(rstandard(model)$pvalue1 0.602$observed.runs1 18$expected.runs1 16.09677$n11 13$n21 18$k 1 0P 值为 0.602，不拒绝误差项的独立性。(d) win.graph(width=3,height=3,pointsize=8)qqnorm(residuals(model)残差的正态 QQ 图selitnauQelpmaSNormal Q-Q Plot-1 0 1Theoretical Quantiles2QQ 图看出有一点小的曲率，但是这种现象可能是俩个极端值造成的。(e)对残差的正态性进行 shapir

14、o-wilk 检验shapiro.test(residuals(model)Shapiro-Wilk normality testdata: residuals(model)W = 0.93509, p-value = 0.06043结果表明，我们不会拒绝通常意义水平的正态性。第九章 预测9.2(a) Y 2007 (1) 5 1.1Y2007 0.5Y2006 5 1.1(10) 0.5(11) 10.5所以， Y 2007(2)5 1.1Y2008 0.5Y2007 5 1.1(10.5) 0.5(10) 11.55(b)从上式中看出 1 1.1 ，1 - 1 0 0, 0 1, 11.1

15、(c) Y2007 (1) 2 210.5 2 2 10.5 2.83 即 2008 年 预测的 95%预测 极限 为7.67-13.33.(d)因为有， Yt 1 (t) Yt1(t 1) t Yt 1 Yt(1)所以， Y 2008(1) Y 2007 (2) 1 Y2008 Y 2007(1) 11.55 1.112 10.5 13.2第十章 季节模型10.12(a) data(boardings);series=boardings,1 plot(series,type='l',ylab='Light Rail&Bus Boardings') po

16、ints(series,x=time(series),pch=as.vector(season(series)roda.21suB.55uB.52liB&1itaR.05th21ihLg1ONFMJ(b) acf(as.vector(series),ci.type='ma')JFAMAONFMNJFMMAA JJ D2001 2002 2003 2004 2005 2006TimeSeries as.vector(series)Lag滞后期为 1，5， 6， 12 ，时候，存在显著的自相关。(c) model=arima(series,order=c(0,0,3),se

17、asonal=list(order=c(1,0,0),period=12);modelCall: arima(x = series, order = c(0, 0, 3), seasonal = list(order = c(1, 0, 0), period = 12)Coefficients:ma1ma2ma3sar1intercept0.72900.61160.29500.877612.5455s.e. 0.11860.11720.11180.05070.0354sigma2 estimated as 0.0006542: log likelihood = 143.54, aic = -2

18、77.09 所有的变量都是显著的。(d) model2=arima(series,order=c(0,0,4),seasonal=list(order=c(1,0,0),period=12);model2 Call:arima(x = series, order = c(0, 0, 4), seasonal = list(order = c(1, 0, 0), period = 12) Coefficients:ma1ma2ma3ma4sar1intercept0.72770.66860.42440.14140.891812.5459s.e. 0.12120.13270.16810.12280

19、.04450.0419sigma2 estimated as 0.0006279: log likelihood = 144.22, aic = -276.45模型 2 中 AIC=-276.45，模型 1 中的 AIC= -277.09，AIC越小越好，所以模型是过度拟合的。 第十二章 异方差时间序列模型12.1 library(TSA)data(CREF)r.cref=diff(log(CREF)*100win.graph(width = 4.875,height = 2.5,pointsize = 8)plot(abs(r.cref)win.graph(width = 4.875,hei

20、ght = 2.5,pointsize = 8)plot(r.cref2)12.9(a) > data(google)> plot(google)收益率数据的时间序列图> acf(google)> pacf(google)根据 ACF和 PACF可以得知，无自相关。(b) 计算 google 日收益率均值。> t.test(google,alternative = 'greater')One Sample t-testdata: googlet = 2.5689, df = 520, p -value = 0.00524alternative hy

21、pothesis: true mean is greater than 095 percent confidence interval:0.000962967 Inf sample estimates:mean of x 0.002685589x的均值为 0.002685589 ，备择假设为： 均值异于 0，根据 P值显示， 0.00524 接受备择假设。(c)McLeod-Li 检验 ARCH效应。> win.graph(width = 4.875,height = 3,pointsize = 8) > McLeod.Li.test(y=google)根据图显示，所有滞后值在 5

22、%的水平上均显著。说明数据具有ARCH特征。(d) 识别 GARCH模型，估计识别的模型并对拟合的模型进行模型诊断检验。> eacf(google2) 取值平方的样本 EACF AR/MA0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 x x o o o o o o o x o o o x 1 x o o o o o o o o x o o o x 2 x o o o o o o o o x o o o x 3 x x x o o o o o o x o o o x 4 x x x o o o o o o o o o o o 5 x x x o o o o o o o

23、 o o o o 6 x x x x o o o o o o o o o o 7 o x x o o x o o o o o o o o> eacf(abs(google) 绝对值的样本 EACF AR/MA0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 x x x o o o x o o x o o x x1 x o o o o o o o o o o o o x2 x x o o o o o o o o o o o x3 x x x o o o o o o o o o o x4 x o x o o o o o o o o o o o5 x o x o x o o

24、o o o o o o o6 o x x x x x o o o o o o o o7 x o x x x o x o o o o o o o根据上图， 得知设定 GARCH(1,1) 模型。 Google 日收益率的平方值相应的样本 EACF 也得知模 型 GARCH(1,1) 符合。检验模型：> m1=garch(x=google,order = c(1,1)> summary(m1)Call:garch(x = google, order = c(1, 1)Model:GARCH(1,1)Residuals:Min 1Q Median 3Q Max-3.64597 -0.46

25、486 0.08232 0.65379 5.73937Coefficient(s):Estimate Std. Error t value Pr(>|t|)a0 5.058e-051.232e-054.106 4.03e-05 *a1 1.264e-012.136e-025.920 3.21e-09 *b1 7.865e-013.578e-0221.980 < 2e-16 *Signif. codes:0 * ' 0.001 * ' 0.01* ' 0.05 . ' 0.1' 1Diagnostic Tests:Jarque Bera Testdata: ResidualsX-squared = 223.86, df = 2, p-value < 2.2e-16Box-Lj