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

6、
方法一：趋势拟合法
income<-scan('习题4.6数据.txt')
ts.plot(income)
由时序图可以看出，该序列呈现二次曲线的形状。

于是，我们对该序列进行二次曲线拟合：
t<-1:length(income)
t2<-t^2
z<-lm(income~t+t2)
summary(z)
lines(z$fitted.values, col=2)
方法二：移动平滑法拟合
选取N=5
income.fil<-filter(income,rep(1/5,5),sides=1)
lines(income.fil,col=3)
7、（1）
milk<-scan('习题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(z.1)
lines(z.s$trend,col=3)
z.2<-ts(z.1)
t<-1:length(z.2)
t2<-t^2
t3<-t^3
r1<-lm(z.2~t)
r2<-lm(z.2~t+t2)
r3<-lm(z.2~t+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<-pt^2
pt3<-pt^3
pt<-matrix(c(pt1,pt2,pt3),byrow=T,nrow=3)/*为预测时间的矩阵。

*/ p<-r3$coef[2:4]%*%pt+r3$coef[1]/*矩阵的乘法为%*%；coef【1】为其截距项，coef【2：4】为其系数*/
p1<-z.s$sea[1: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<-diff(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）模型，故综合可采用乘积模型12(0,1,1)(0,1,1)SARIMA ##即ri1、ma1模型乘以季节因素
result<-arima(z,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))/*季节因素里的order 为阶数的意思，与前面的airma 模型的阶数含义同*/
tsdiag(result)//诊断
##下图为预测后的图
4.8
z<-scan('4.8.txt')
adf.test(z) ##单位根检验。

比较科学的定量的方法
##其原假设：具有单位根，即不平稳。

此题中接受备则假设：平稳。

指数平滑预测
ffe<-function(z,a) ##定义指数平滑预测。

其中a为平滑项
{
y<-c()
y<-z[1]
for(i in 1:length(z))
y<-c(y, a*z[i]+(1-a)*y[i])
return(y)
}
y<-ffe(z,0.6) ##执行上述定义的function
ts.plot(z)
lines(y,col=3)
y[length(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.1[length(z.1)]
z.2<-z[(length(z)-10):length(z)]
for(i in 1:n)
{
m<-sum(rep(1/12,12-i)*z.2[i: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)
##SARIMA
par(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),period=12
),fixed=c(NA,NA,0,NA,NA,NA)) ##观察图，发现第三项在置信区间内，故认为可能为限定的sarima模型。

最后两个NA指季节指数中的sar1和sma1.
##第三个的aic值最小，即模型拟合效果最好
tsdiag(result) ##检验通过。