kings<-scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)
#读入时间序列数据,忽略前三行
kings
kingstimeseries<-ts(kings)#将数据存储到一个时间序列对象中去
births<-scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")
births
birthstimeseries<-ts(births,frequency=12,start=c(1946,1))
birthstimeseries
souvenir<-scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries<-ts(souvenir,frequency=12,start=c(1987,1))
souvenirtimeseries
plot.ts(kingstimeseries)#可用相加模型来描述,
plot.ts(birthstimeseries)#
plot.ts(souvenirtimeseries)#相加模型不适合描述,季节性波动和随机变动的大小随时间序列逐步上升
logsouvenirtimeseries<-log(souvenirtimeseries)#自然对数转换
plot.ts(logsouvenirtimeseries)
library("TTR")
kingstimeseriesSMA3<-SMA(kingstimeseries,n=3)#跨度为3的简单移动平均平滑
plot.ts(kingstimeseriesSMA3)
kingstimeseriesSMA8<-SMA(kingstimeseries,n=8)#跨度为3的简单移动平均平滑
plot.ts(kingstimeseriesSMA8)
birthstimeseriescomponents<-decompose(birthstimeseries)
birthstimeseriescomponents$seasonal#季节
birthstimeseriescomponents$trend#趋势
birthstimeseriescomponents$random#随机
plot(birthstimeseriescomponents)
birthstimeseriesseasonallyadjusted<-birthstimeseries-birthstimeseriescomponents$seasonal
plot(birthstimeseriesseasonallyadjusted)
rain<-scan("http://robjhyndman.com/tsdldata/hurst/precip1.dat",skip=1)
rainseries<-ts(rain,start=c(1813))
plot.ts(rainseries)
rainseriesforecasts<-HoltWinters(rainseries,beta=FALSE,gamma=FALSE)
rainseriesforecasts
rainseriesforecasts$fitted
plot(rainseriesforecasts)
rainseriesforecasts$SSE
library("forecast")
rainseriesforecasts2<-forecast.HoltWinters(rainseriesforecasts,h=8)
plot.forecast(rainseriesforecasts2)
rainseriesforecasts2$residuals
acf(rainseriesforecasts2$residuals,lag.max=20)
plotForecastErrors<-function(forecasterrors)
{
mybinsize<-IQR(forecasterrors)/4
mysd<-sd(forecasterrors)
mymin<-min(forecasterrors)-mysd*3
mymax<-max(forecasterrors)+mysd*3
mybins<-seq(mymin,mymax,mybinsize)
hist(forecasterrors,col="red",freq=FALSE,breaks=mybins)
mynorm<-rnorm(10000,mean=0,sd=mysd)
myhist<-hist(mynorm,plot=FALSE,breaks=mybins)
points(myhist$mids,myhist$density,type="l",col="blue",lwd=2)
}
plotForecastErrors(rainseriesforecasts2$residuals)
霍尔特指数平滑法
skirts<-scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
skirtsseries<-ts(skirts,start=c(1866))
plot.ts(skirtsseries)
skirtsseriesforecasts<-HoltWinters(skirtsseries,gamma=FALSE)
skirtsseriesforecasts
skirtsseriesforecasts$SSE
plot(skirtsseriesforecasts)
skirtsseriesforecasts2<-forecast.HoltWinters(skirtsseriesforecasts,h=19)
plot.forecast(skirtsseriesforecasts2)
acf(skirtsseriesforecasts2$residuals,lag.max=20)
Box.test(skirtsseriesforecasts2$residuals,lag=20,type="Ljung-Box")
plot.ts(skirtsseriesforecasts2$residuals)
plotForecastErrors(skirtsseriesforecasts2$residuals)
Holt-Winters指数平滑法
souvenir<-scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries<-ts(souvenir,frequency=12,start=c(1987,1))
souvenirtimeseries
logsouvenirtimeseries<-log(souvenirtimeseries)
souvenirtimeseriesforecasts<-HoltWinters(logsouvenirtimeseries)
souvenirtimeseriesforecasts
souvenirtimeseriesforecasts$SSE
plot(souvenirtimeseriesforecasts)
souvenirtimeseriesforecasts2<-forecast.HoltWinters(souvenirtimeseriesforecasts,h=48)
plot.forecast(souvenirtimeseriesforecasts2)
acf(souvenirtimeseriesforecasts2$residuals,lag.max=20)
Box.test(souvenirtimeseriesforecasts2$residuals,lag=20,type="Ljung-Box")
plot.ts(souvenirtimeseriesforecasts2$residuals)
plotForecastErrors(souvenirtimeseriesforecasts2$residuals)
###################################################################
###ARMA 指数平滑法 ###
#####################################################################
skirtsseriesdiff1<-diff(skirtsseries,differences=1)
plot.ts(skirtsseriesdiff1)
skirtsseriesdiff2<-diff(skirtsseries,differences=2)
plot.ts(skirtsseriesdiff2)
kings<-scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)#读入时间序列数据,忽略前三行
kings
kingstimeseries<-ts(kings)#将数据存储到一个时间序列对象中去
kingstimeseriesdiff1<-diff(kingstimeseries,differences=1)
plot.ts(kingstimeseriesdiff1)
acf(kingstimeseriesdiff1,lag.max=20)#画出相关图
acf(kingstimeseriesdiff1,lag.max=20,plot=FALSE)#取得相关系数
pacf(kingstimeseriesdiff1,lag.max=20)#画出相关图
pacf(kingstimeseriesdiff1,lag.max=20,plot=FALSE)#取得相关系数
auto.arima(kings)
kingstimeseriesarima<-arima(kingstimeseries,order=c(0,1,1))
kingstimeseriesarima
kingstimeseriesforecasts<-forecast.Arima(kingstimeseriesarima,h=5,level=c(99.5))
kingstimeseriesforecasts
plot.forecast(kingstimeseriesforecasts)
acf(kingstimeseriesforecasts$residuals,lag.max=20)
Box.test(kingstimeseriesforecasts$residuals,lag=20,type="Ljung-Box")
plot.ts(kingstimeseriesforecasts$residual)
plotForecastErrors(kingstimeseriesforecasts$residual)
volcanodust<-scan("http://robjhyndman.com/tsdldata/annual/dvi.dat",skip=1)
volcanodustseries<-ts(volcanodust,start=c(1500))
plot.ts(volcanodustseries)
acf(volcanodustseries,lag.max=20)#画出相关图
acf(volcanodustseries,lag.max=20,plot=FALSE)#取得相关系数
pacf(volcanodustseries,lag.max=20)#画出相关图
pacf(volcanodustseries,lag.max=20,plot=FALSE)#取得相关系数
auto.arima(volcanodustseries)
volcanodustseriesarima<-arima(volcanodustseries,order=c(2,0,0))
volcanodustseriesarima
volcanodustseriesforecasts<-forecast.Arima(volcanodustseriesarima,h=31)
volcanodustseriesforecasts
plot.forecast(volcanodustseriesforecasts)
acf(volcanodustseriesforecasts$residuals,lag.max=20)
Box.test(volcanodustseriesforecasts$residuals,lag=20,type="Ljung-Box")
plot.ts(volcanodustseriesforecasts$residual)
plotForecastErrors(volcanodustseriesforecasts$residual)
mean(volcanodustseriesforecasts$residuals)
#