2018.10.09R-SARIMA学习:程序

2018.10.09R-SARIMA学习

library(forecast)
library(tseries)
data(AirPassengers)
class(AirPassengers)
start(AirPassengers)
end(AirPassengers)
frequency(AirPassengers)
summary(AirPassengers)
tsair <- ts(AirPassengers, start=c(1949, 1), frequency=12)
tsair
tsair.subset<- window(tsair,start=c(1949,1),end=c(1960,6),frequency=12)
tsair.subset
plot(AirPassengers,ylab="AirPassengers")
plot(tsair.subset,ylab="AirPassengers")
plot(log(AirPassengers),ylab="log(AirPassengers)")
stl(tsair.subset, s.window="periodic", t.window=)#得到是取对数后的序列分解趋势季节随机组成
fit1 <- stl(tsair.subset, s.window="periodic", t.window=)
exp(fit1$time.series)#exp(fit$time.series)回到原序列
plot(fit1)
#--------------建模
ndiffs(log(tsair.subset))
ndiffs(log(AirPassengers))
dla<-diff(log(tsair.subset),d=1,D=1)
dla2<-diff(log(AirPassengers),d=1) 
plot(dla)
plot(dla2)
adf.test(dla)#alternative hypothesis:stationary, p-value=0.01<0.05,stationary
Acf(dla) #q=1,Q=0,1
Pacf(dla)#p=1,2,P=0,1

D <- arima(dla,c(1,1,1), seasonal = list(order=c(0,1,1),period=12),method = "ML",include.mean = FALSE)
D1 <- arima(dla2, c(1,1,1), seasonal= list(order=c(0,1,1),period=12),method = "ML",include.mean = FALSE)
D       #aic=-442.15  （负数）最小：SARIMA（1,1,1）（0,1,1）12
accuracy(D) #RMSE=0.03584
summary(D)
qqnorm(D$residuals)   
qqline(D$residuals) #拟合结果较好
Box.test(D$residuals, type="Ljung-Box")#p-value=0.72,为白噪声
Box.test(D$residuals, type="Ljung-Box",lag=6)#p-value=0.2362
Box.test(D$residuals, type="Ljung-Box",lag=12)#p-value=0.4703
Box.test(D$residuals, type="Ljung-Box",lag=1)#p-value=0.7665
pre=predict(D,n.ahead = 6)
U=pre$pred+1.96*pre$se
L=pre$pred-1.96*pre$se
ts.plot(pre$pred,xlab = 'Year',
        ylab = 'log(Air Passengers)',
        main="ARIMA(1,1,1)(0,1,1)12预测的1960年7月-12月的飞机乘客对数量",col="red")
lines(U,col="blue")
lines(L,col="blue")
forecast(D, 6)
a <- forecast(D,7)
b = forecast(D,1)
plot(a,xlab = 'Year',
     ylab = 'log(Air Passengers)',
     main="ARIMA(1,1,1)(0,1,1)12预测的1960年7月-12月的飞机乘客对数量",col="red")
par(new=TRUE)
plot(b,xlab = 'Year',
     ylab = 'log(Air Passengers)',
     main="ARIMA(1,1,1)(0,1,1)12预测的1960年7月-12月的飞机乘客对数量")

;

LUFEI第三周学习笔记

1.SARIMA(Seasonal auto regression integrated moving average)模型
一些时间序列存在季节性周期波动，需加入季节性算子，被称为SARIMA模型（或季节ARIMA模型），SARIMA模型是随机季节模型与ARIMA模型的结合，表示为SARIMA(p,d,q)(P,D,Q)s。

1.1建模步骤
①序列平稳化：非季节差分与季节差分
②模型识别：确定SARIMA(p,d,q)(P,D,Q)s模型中的参数p,d,q,P,D,Q,s的值。
③模型的参数估计和模型诊断：参数t检验有无统计学意义；AIC/BIC值选择模型
④模型预测：以预测值与实际值的相对误差评价模型的预测效果。

1.2非平稳时间序列
时间序列分为平稳性的与非平稳性的。非平稳时间序列可能包含四种趋势：长期趋势(T)、循环趋势©、季节趋势(S)和不规则趋势(I)。
①长期趋势T：是时间序列在长时期内呈现出来的持续向上或持续向下的变动。
②季节变动S：是时间序列在一年内重复出现的周期性波动。
③循环波动C：是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间，但与趋势不同，它不是朝着单一方向的持续变动，而是涨落相同的交替波动。
④不规则波动I：是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中，致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列

1.2.1时间序列季节调整方法主要有四种趋势模型：加法模型、乘法模型、对数加法模型和伪加法模型。
①加法模型主要适用于呈线性增长的数据序列，或者是围绕某一个中值波动的数据序列
②乘法模型主要适用于呈指数级数增长的序列
③对数加法模型主要适用于同比增速呈线性增长的数据序列
④伪加法模型则主要是对某些非负时间序列进行季节调整，它们具有这样的性质：在每一年中的相同月份出现接近于0的正值，在这些月份含有接近于0的季节因子，受这些小因子的影响，季节调整结果将出现偏差，如：农产品产量

1.3常用的模型有乘法模型和加法模型
1.3.1乘法模型和加法模型都是将形成时间序列变动的四类构成因素，按照它们的影响方式不同，设定的组合模型：
乘法模型: Y = T·S·C·I 加法模型: Y = T+S+C+I
式中:Y:时间序列的指标数值 T:长期趋势 S:季节趋势 C:循环趋势 I:不规则趋势

1.3.2两个模型季节因素的表述的不同在于:
①乘法模型是假定四个因素对现象的发展的影响是相互作用的，用于相对数总变量的计算。以长期趋势成分的绝对量为基础，其余量均以比率表示。
②加法模型是假定四个因素的影响是相互独立的，用于总量指标总变动的计算。每个成分均以绝对量表示。
假设我们有一个时序，记录了10年来摩托车的月销量。在可加模型中， 11月和12月（圣诞节）的销量一般会增加500，而1月（一般是销售淡季）的销量则会减少200。此时季节性的波动量和当时的销量无关。在相乘模型中，11月和12月的销量则会增加20%， 1月的销量减少10%，即季节性的波动量和当时的销量是成比例的。很多时候，相乘模型比相加模型更现实一些。

1.4在R中实现SARIMA建模
library(forecast)
library(tseries)
data(AirPassengers)
class(AirPassengers)
start(AirPassengers)
end(AirPassengers)
frequency(AirPassengers)
summary(AirPassengers)

tsair <- ts(AirPassengers, start=c(1949, 1), frequency=12)
Tsair

tsair.subset<- window(tsair,start=c(1949,1),end=c(1960,6),frequency=12)
tsair.subset

plot(AirPassengers,ylab=“AirPassengers”)

plot(tsair.subset,ylab=“AirPassengers”)

plot(log(AirPassengers),ylab=“log(AirPassengers)”)

stl(tsair.subset, s.window=“periodic”, t.window=)#得到是取对数后的序列分解趋势季节随机组成
fit1 <- stl(tsair.subset, s.window=“periodic”, t.window=)
exp(fit1KaTeX parse error: Expected 'EOF', got '#' at position 13: time.series)#̲exp(fittime.series)回到原序列

plot(fit1)

#--------------建模
ndiffs(log(tsair.subset))
ndiffs(log(AirPassengers))

dla<-diff(log(tsair.subset),d=1,D=1)
dla2<-diff(log(AirPassengers),d=1)
plot(dla)

plot(dla2)

adf.test(dla)#alternative hypothesis:stationary, p-value=0.01<0.05,stationary

Acf(dla) #q=1,Q=0,1

Pacf(dla)#p=1,2,P=0,1

fit2 <- arima(dla,c(1,1,1), seasonal = list(order=c(0,1,1),period=12))
fit0 <- arima(dla2, c(1,1,1), seasonal= list(order=c(0,1,1),period=12))
fit2 #aic=-442.15 （负数）最小：SARIMA（1,1,1）（0,1,1）12

accuracy(fit2) #RMSE=0.03584

fit3 <- arima(dla,c(2,1,1), seasonal = list(order=c(0,1,0),period=12))
fit3 #aic=-411.59
accuracy(fit3) #RMSE=0.041465
fit4 <- arima(dla,c(1,1,1), seasonal = list(order=c(0,1,0),period=12))
fit4 #aic=-413.58
accuracy(fit4) #RMSE=0.041463
fit5 <- arima(dla,c(2,1,1), seasonal = list(order=c(0,1,1),period=12))
fit5 #aic=-440.4
accuracy(fit5) #RMSE=0.035793
fit6 <- arima(dla,c(1,1,1), seasonal = list(order=c(1,1,1),period=12))
fit6 #aic=-440.35
accuracy(fit6) #RMSE=0.03581
fit7 <- arima(dla,c(2,1,1), seasonal = list(order=c(1,1,0),period=12))
fit7 #aic=-435.3
accuracy(fit7) #RMSE=0.03681
fit8 <- arima(dla,c(1,1,1), seasonal = list(order=c(1,1,0),period=12))
fit8 #aic=-436.84
accuracy(fit8) #RMSE=0.03692
fit9 <- arima(dla,c(2,1,1), seasonal = list(order=c(1,1,1),period=12))
fit9 #aic=-438.72
accuracy(fit9) #RMSE=0.03574

qqnorm(fit2$residuals)

qqline(fit2$residuals) #拟合结果较好

Box.test(fit2$residuals, type=“Ljung-Box”) #p-value=0.72,为白噪声

qqnorm(fit5$residuals) qqline(fit5$residuals)qqline(fit5residuals)
forecast(fit2, 6)

a <- forecast(fit2,7)
plot(a,xlab = ‘Year’,
ylab = ‘log(Air Passengers)’,
main=“ARIMA(1,1,1)(0,1,1)12预测的1960年7月-12月的飞机乘客对数量”)

par(new=TRUE)
plot(b,xlab = ‘Year’,
ylab = ‘log(Air Passengers)’,
main=“ARIMA(1,1,1)(0,1,1)12预测的1960年7月-12月的飞机乘客对数量”)