R语言时间序列作业.docx

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 表示只显示各数据点之间的连接线段；type=p 只想显示数据点。points(y=hours,x=time(hours),pch=as.vector(season(hours)3.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: Min 1Q Median 3Q Max -1.00603 -0.25431 -0.02267 0.22884 0.98358 Coefficients: 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 * 0.01 * 0.05 . 0.1 1Residual standard error: 0.423 on 57 degrees of freedomMultiple R-squared: 0.5921, Adjusted R-squared: 0.5778 F-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(rstudent(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(hours.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 第五章 非平稳时间序列模型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.5 2=-0.5 1=0.5 2=-0.255.7(a) A:AR(2) 1=0.9 2=0.09 B: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=Winnebago 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=Percentage Change,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)(b) acf(series)第八章 模型诊断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(rstandard(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: residuals(model)W = 0.93509, p-value = 0.06043结果表明，我们不会拒绝通常意义水平的正态性。第九章 预测9.2(a)所以，(b)从上式中看出，(c)即2008年预测的95%预测极限为7.67-13.33.(d)因为有，所以，第十章 季节模型10.12(a) data(boardings);series=boardings,1 plot(series,type=l,ylab=Light Rail&Bus Boardings) points(series,x=time(series),pch=as.vector(season(series)(b) acf(as.vector(series),ci.type=ma)滞后期为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: ma1 ma2 ma3 sar1 intercept 0.7290 0.6116 0.2950 0.8776 12.5455s.e. 0.1186 0.1172 0.1118 0.0507 0.0354sigma2 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: ma1 ma2 ma3 ma4 sar1 intercept 0.7277 0.6686 0.4244 0.1414 0.8918 12.5459 s.e. 0.1212 0.1327 0.1681 0.1228 0.0445 0.0419sigma2 estimated as 0.0006279: log likelihood = 144.22, aic = -276.45模型2中AIC=-276.45，模型1中的AIC= -277.09，AIC越小越好，所以模型是过度拟合的。第12章 异方差时间序列模型12.1 library(TSA) data(CREF) r.cref=diff(log(CREF)*100 win.graph(width = 4.875,height = 2.5,pointsize = 8) plot(abs(r.cref) win.graph(width = 4.875,height = 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 hypothesis: true mean is greater than 095 percent confidence interval: 0.000962967 Infsample estimates: mean of x 0.002685589 x的均值为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%的水平上均显著。说明数据具有ARCH特征。(d)识别GARCH模型，估计识别的模型并对拟合的模型进行模型诊断检验。 eacf(google2)取值平方的样本EACFAR/MA 0 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 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 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)绝对值的样本EACFAR/MA 0 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 x 1 x o o o o o o o o o o o o x 2 x x o o o o o o o o o o o x 3 x x x o o o o o o o o o o x 4 x o x o o o o o o o o o o o 5 x o x o x o o o o o o o o o 6 o x x x x x o o o o o o o o 7 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.46486 0.08232 0.65379 5.73937 Coefficient(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 Be