时间序列分析王燕习题4答案2

1、6、方法一：趋势拟合法income<-scan('习题4.6数据.txt')ts.plot(income)由时序图可以看出，该序列呈现二次曲线的形状。于是，我们对该序列进行二次曲线拟合：t<-1:length(income)t2<-t2z<-lm(incomet+t2)summary(z)lines(z$fitted.values, col=2)方法二：移动平滑法拟合选取N=5income.fil<-filter(income,rep(1/5,5),sides=1)lines(income.fil,col=3)7、（1）milk<-scan(

2、'习题4.7数据.txt')ts.plot(milk)从该序列的时序图中，我们看到长期递增趋势和以年为固定周期的季节波动同时作用于该序列,因此我们可以采用乘积模型和加法模型。在这里以加法模型为例。z<-scan('4.7.txt')ts.plot(z)z<-ts(z,start=c(1962,1),frequency=12)z.s<-decompose(z,type='additive') /运用加法模型进行分解z.1<-z-z.s$seas /提取其中的季节系数，并在z中减去（因为是加法模/型）该季节系数ts.plot(

3、z.1)lines(z.s$trend,col=3)z.2<-ts(z.1)t<-1:length(z.2)t2<-t2t3<-t3r1<-lm(z.2t)r2<-lm(z.2t+t2)r3<-lm(z.2t+t2+t3)summary(r1)summary(r2)summary(r3) #发现3次拟合效果最佳，故选用三次拟合ts.plot(z.2)lines(r3$fitt,col=4)pt<-(length(z.2)+1) : (length(z.2)+12)pt1<-pt #预测下一年序列pt2<-pt2pt3<-pt3p

4、t<-matrix(c(pt1,pt2,pt3),byrow=T,nrow=3)/*为预测时间的矩阵。*/p<-r3$coef2:4%*%pt+r3$coef1/*矩阵的乘法为%*%；coef【1】为其截距项，coef【2：4】为其系数*/p1<-z.s$sea1:12+p/*加回原有季节系数，因为原来是加法模型*/ts.plot(ts(z),xlim=c(1,123),ylim=c(550,950)lines(pt1,p1,col=2)#包含季节效应的 SARIMA模型z<-scan('4.7.txt')ts.plot(diff(z)sq<-di

5、ff(diff(z),lag=12) /*12步差分*/par(mfrow=c(2,1) acf(sq,50) pacf(sq,50)#观察上图，发现ACF图12阶处明显，24阶处即变到置信区间内。#而PACF图12阶，24阶，36阶处有一个逐渐递减过程，可认为#拖尾，故可以考虑对季节效应部分采用MA(1)模型#同时，ACF图在第一阶处显著后即立刻变动到置信区间内，具有#截尾性质，PACF图在第5、6阶时变动到置信区间外，可以考虑#使用MA（1）模型，故综合可采用乘积模型#即ri1、ma1模型乘以季节因素result<-arima(z,order=c(0,1,1),seasonal=li

6、st(order=c(0,1,1),period=12)/*季节因素里的order为阶数的意思，与前面的airma模型的阶数含义同*/tsdiag(result)/诊断#下图为预测后的图4.8z<-scan('4.8.txt')adf.test(z) #单位根检验。比较科学的定量的方法#其原假设：具有单位根，即不平稳。此题中接受备则假设：平稳。指数平滑预测ffe<-function(z,a) #定义指数平滑预测。其中a为平滑项y<-c()y<-z1for(i in 1:length(z)y<-c(y, a*zi+(1-a)*yi)return(y)

7、y<-ffe(z,0.6) #执行上述定义的function ts.plot(z) lines(y,col=3) ylength(y)简单移动平均z.1<-filter(z,rep(1/12,12),side=1) #side=1是指将所有算不出的序列值都空到最前面去，而在尾部没有空值。z.1<-c(NA,z.1)ts.plot(z)lines(z.1,col=3)meand<-function(z,z.1,n) #预测函数。以12为周期。依次为原始数据，平滑值，预测步数y<-z.1length(z.1)z.2<-z(length(z)-10):length

8、(z)for(i in 1:n) m<-sum(rep(1/12,12-i)*z.2i:length(z.2) n<-sum(rep(1/12,i)*y) y<-c(y,m+n) #一直重复：预测，原始数列取代一个，预测数列拿来一个return(y)y<-meand(z,z.1,11)y<-c(z.1,y) ts.plot(z,xlim=c(0,205) lines(y,col=3)#SARIMApar(mfrow=c(2,1)ds<-diff(z)acf(ds,40)pacf(ds,40)#可以看出有一些不明显的周期性，故采用sarima拟合result<-arima(z,order=c(2,1,0),seasonal=list(order=c(1,0,0),period=12)#在季节部分很少出现2以上的数字（指seasonal中的order部分） result<-arima(z,order=c(2,1,0),seasonal=list(order=c(1,0,1),period=12) result<-arima(z,order=c(4,1,0),seasonal=list(order=c(1,0,1