R语言时间序列作业.docx

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

2、s)CT UOH wrkvnoM7 09605094093 03.10(a) data(hours);hours.lm=lm(hourstime(hours)+I(time(hours)A2);summary(hours.lm)用最小二乘法拟合二次趋势，结果显示如下：Call: lm(formula = hours time(hours) + I(time(hours)A2)Residuals:Min 1Q Median 3Q Max -1.00603 -0.25431 -0.02267 0.22884 0.98358Coefficients:Estimate Std. Error t val

3、ue 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)A2) -1.299e-01 2.933e-02 -4.428 4.35e-05 * Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.423 on 57 degrees of freedom Multiple R-squared: 0.5921, Adju

4、sted 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 Residuals) points(y=rstudent(hours.lm),x=as.vector(time(hours),pch=as.vector(season(hours)标准残差的时间序列，应用月度绘图标志。(为了更容易识别季节性)带季节性图标的的残差-时间图(c)runs(rstud

5、ent(hours.lm)对标准差进行游程检验$pvalue1 0.00012$observed.runs1 16$expected.runs1 30.96667$n11 31$n21 29$k1 0结果解释：P值为0.00012,表明非随机性是合理的。(d)acf(rstudent(hours.lm)标准残差的样本自相关函数(e) qqnorm(rstudent(hours.lm);qqline(rstudent(hours.lm)(QQ 图)正态性可以通过正态得分或者分位数 -分位数(QQ)图来检验。此处的直线型图形支持了该 模型中随机项是正态分布的假设。hIst(rstudent(hou

6、rs.lm),xlab=StandardIzed Residuals)标准残差的直方图(季节均值模型的标准残差直方图)shapiro.test(rstudent(hours.lm)正态性检验(Shapiro-Wilk检验)本质是：计算残差与相应的正态分位数之间的相关系数。相关性越小，就越有理由否定正态性。Shapiro-Wilk normality testdata: rstudent(hours.lm)W = 0.99385, p-value = 0.9909根据上面的检验结果，我们不能拒绝模型的随机项是正态分布的假设。第四章平稳时间序列模型4.4 1 1第五章非平稳时间序列模型5.1(a)

7、 ARMA (2,1) p=2,q=1,参数值 4 和。4 1=1 力 2=-0.25 1=0.1(b)IMA(2,0) p=2,d=1,q=0,参数值 4 和 0(c)ARMA(2,2) p=2,q=2,参数值 4 和 0()1=0.5 ()2=-0.5 1=0.5 2=-0.255.7(a) A:AR(2) 力 1=0.9 4 2=0.09B:IMA(1,1) 1=0.1(b)一个是固定的一个是不固定的。5.11(a) data(winnebago);win.graph(width=6.5,height=3,pointsize=8)plot(winnebago,type=o,ylab=Wi

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

9、hange,xlab=Difference of Logs)cor(diff(log(winnebago)-1,percentage-1)1 0.9646886结果显示：0.96认为一致。第六章模型识别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)Lag(

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

11、d(model)$pvalue1 0.602$observed.runs1 18$expected.runs1 16.09677$n11 13$n21 18$k1 0P值为0.602,不拒绝误差项的独立性。(d)win.graph(width=3,height=3,pointsize=8)qqnorm(residuals(model)残差的正态QQ图QQ图看出有一点小的曲率，但是这种现象可能是俩个极端值造成的。(e)对残差的正态性进行 shapiro-wilk检验shapiro.test(residuals(model)Shapiro-Wilk normality testdata: resi

12、duals(model)W = 0.93509, p-value = 0.06043结果表明，我们不会拒绝通常意义水平的正态性。第九章预测所以，Y 2008 (1) Y 2007 (2)9.2(a)Y2007(1) 5 1.1Y2007 0.5Y2006 5 1.1(10) 0.5(1 1) 10.5所以，Y 2007 (2)51 .1Y20080.5Y2007 51.1(10.5) 0.5(10) 11.55(b)从上式中看出11.1 ,1-10 0,0 1, 1 1 1.1(c) Y2007 (1) 2 .210.52.2 10.52.83即2008年预测的95%预测极限为7.67-13.

13、33.(d)因为有,Yt 1(t)丫奇1) t Y 1丫1 Y2008 Y2007 (1) 11.55 1.1 12 10.5 13.2第十章季节模型10.12(a) data(boardings);series=boardings,1plot(series,type=T,ylab=Light Rail&Bus Boardings)points(series,x=time(series),pch=as.vector(season(series)aB sB2001 2002 2003 2004 2005 2006Time(b) acf(as.vector(series),ci.type=ma)滞

14、后期为1 , 5, 6, 12,时候，存在显著的自相关。(c) model=arima(series,order=c(0,0,3),seasonal=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.0

15、354sigmaA2 estimated as 0.0006542: log likelihood =143.54, aic = -277.09所有的变量都是显著的。(d)model2=arima(series,order=c(0,0,4),seasonal=list(order=c(1,0,0),period=12);model2Call:arima(x = series, order = c(0, 0, 4), seasonal = list(order = c(1,0, 0), period = 12) Coefficients:ma1ma2ma3ma4sar1intercept0.72

16、770.66860.42440.14140.891812.5459s.e.0.12120.13270.16810.12280.04450.0419sigmaA2 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,p

17、ointsize = 8)plot(abs(r.cref)win.graph(width = 4.875,height = 2.5,pointsize = 8) plot(r.crefA2)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

18、.00524alternative hypothesis: true mean is greater than 095 percent confidence interval:0.000962967 Infsample estimates:mean of x0.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)以-u 心 “ - e 廿

19、弋- u m : 7 。区 uI H I I9 5 10 15 25Lag根据图显示，所有滞后值在 5%的水平上均显著。说明数据具有 ARCH特征。(d)识别GARCH莫型，估计识别的模型并对拟合的模型进行模型诊断检验。 eacf(googleA2)取值平方的样本 EACFAR/MA0 1 2 3 4 5 6 7 8 9 10 11 12 130 x x o o o o o o o x o o o x1x o o o o o o o o x o o o x2x o o o o o o o o x o o o x3x x x o o o o o o x o o o x4x x x o o o o

20、 o o o o o o o5x x x o o o o o o o o o o o6x x x x o o o o o o o o o o7o x x o o x o o o o o o o oeacf(abs(google)绝对值的样本 EACFAR/MA0 1 2 3 4 5 6 7 8 9 10 11 12 130 x x x o o o x o o x o o x x1x o o o o o o o o o o o o x2x x o o o o o o o o o o o x3x x x o o o o o o o o o o x4x o x o o o o o o o o o o

21、 o5x o x o x o o o o o o o o o6o x x x x x o o o o o o o o7x 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.46486

22、 0.08232 0.65379 5.73937Coefficient(s):Estimate Std. Error t value Pr(|t|)a0 5.058e-05 1.232e-05 4.106 4.03e-05 *a1 1.264e-01 2.136e-02 5.920 3.21e-09 *b1 7.865e-01 3.578e-02 21.980 2e-16 *Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Diagnostic Tests:Jarque Bera Testdata: ResidualsX-squared = 223.8

23、6, df = 2, p-value plot(residuals(m1),type=h,ylab=Standardized Residuals)标准残差的QQ正态得分图win.graph(width = 2.5,height = 2.5,pointsize = 8)qqnorm(residuals(m1);qqline(residuals(m1)Theoretical Quant ilea如果模型识别正确，那么残差应该是近似独立同分布的。从上 qq图中显示并不明确。所以我们 进行广义混合检验。绘制google 日收益率的 GARCHI (1, 1)模型的标准残差平方的样本 ACFacf(re

24、siduals(m1)A2,na.action=na.omit)Series residuQls(m1A2排SCD0从图形中得出总体印象是残差平方序列不相关。进行广义混合检验得到的 p值gBox(m1,method=squared)结果显示错误！ ！ ！不清楚为什么！ ！ ！在这里继续使用绝对标准残差重新对模型进行检验o绝对标准残差的样本 ACFacf(abs(residuals(m1),na.action=na.omit)gBox(m1,method=absolute)(e)绘制并评论估计的条件方差的时间序列图图中显示出来有几个时期的波动率较高。(f)对拟合模型的标准残差绘制 QQ图qqnorm(residuals(m1);qqline(residuals(m1)残差看起来并不正常。(g)构造 bl 的 95%的置信区间。(0.7865-1.96*0.03578 , 0.7865+1.96*0.03578 ) =(0.7164,0.8566)