时间序列学习记录

#1.用ts()建立时间序列
#例子1
kings<-scan("//robjhyndman.com/tsdldata/misc/kings.dat",skip=3);
kings
kingstimeseries<-ts(kings)
kingstimeseries
plot(kingstimeseries)

#例子2
#使用frequency来设置一些季度、阅读灯不同的间隔或者频数,
#数据集是从 1946 年 1 月到 1959 年 12 月的纽约每月出生人口数量
births <-scan("//robjhyndman.com/tsdldata/data/nybirths.dat")
births
birthstimeseries<-ts(births,start=c(1946,1),frequency=12)
birthstimeseries
#昆士兰海滨度假圣地的纪念品商店从 1987 年 1 月到 1987 年 12 月的每月销售数据
#原始数据源于 Wheelwright and Hyndman，1998）
souvenir <-scan("//robjhyndman.com/tsdldata/data/fancy.dat")
souvenirtimeseries <- ts(souvenir, frequency=12,start=c(1987,1))
souvenirtimeseries


#绘制时间序列使用plot或者plot.ts命令
plot.ts(kings)
plot.ts(kingstimeseries)
plot.ts(birthstimeseries)
plot(birthstimeseries)
plot.ts(souvenirtimeseries)
#为了更好的体现时间序列的关系，有时候需要对元数据进行变换
logsouvenirtimeseries <- log(souvenirtimeseries)
plot.ts(logsouvenirtimeseries)


#2.如何分解时间序列，分解一个时间序列意味着把它拆分成构成元件，一般序列包含一个趋势部分、一个不规则部分，如果是一个季节
#性时间序列，则还有一个季节性部分

#分解非季节性数据
#使用TTR包
install.packages("TTR")
library("TTR")
#使用SMA平滑时间序列
kingstimeseriesSMA3 <-SMA(kingstimeseries,n=3)#每三个值求平均作为新值
plot.ts(kingstimeseriesSMA3)
#上面的图像显示波动仍然比较大，我们增大跨度平滑令n=8
kingstimeseriesSMA8 <- SMA(kingstimeseries,n=8)
kingstimeseries
plot.ts(kingstimeseriesSMA8)
#观察趋势比较明显了，这个过程可能要反复试验多次
plot.ts(kingstimeseriesSMA5 <- SMA(kingstimeseries,n=5))
plot.ts(kingstimeseriesSMA10 <- SMA(kingstimeseries,n=10))


#基本趋势国王的年龄从55降到38，然后升高到73左右


#分解季节性数据
#一个季节性时间序列包含一个趋势部分，一个季节性部分和一个不规则部分。分解时间序列就意味着要把时间序
#列分解称为这三个部分：也就是估计出这三个部分。
#对于可以使用相加模型进行描述的时间序列中的趋势部分和季节性部分，我们可以使用 R 中“decompose()”函
#数来估计。这个函数可以估计出时间序列中趋势的、季节性的和不规则的部分，而此时间序列须是可以用相加模
#型描述的。

#继续使用纽约出生人口数据
#decompose命令
birthstimeseriescomponents <- decompose(birthstimeseries)
#估计出的季节性、趋势的和不规则部分现在被存储在变量birthstimeseriescomponents$seasonal,
#birthstimeseriescomponents$trend 和 birthstimeseriescomponents$random  中
birthstimeseriescomponents$seasonal#估计季节性因素，最大的出现在JUL，低估在FEB
plot(birthstimeseriescomponents)#画出观察趋势图以及分解后各部分的趋势图

#我们也可以使用decompose对季节性因素进行调节，例如去除季节性因素
birthstimeseriesseasonallyadjusted <- birthstimeseries -birthstimeseriescomponents$seasonal#去除季节性因素
plot.ts(birthstimeseriesseasonallyadjusted)



#3.使用指数平滑进行预测，只能进行短期预测
#3.1简单指数平滑法，如果你有一个可用相加模型描述的，并且处于恒定水平和没有季节性变动的时间序列，你可以使用简单指数平滑
#法对其进行短期预测。
#数据是伦敦从 1813 年到 1912 年全部的每年每英尺降雨量（初始数据来自 Hipel and McLeod,1994）
rain<-scan("//robjhyndman.com/tsdldata/hurst/precip1.dat",skip=1)
rainseries<-ts(rain,start=c(1813))
plot.ts(rainseries)

#R中使用HoltWinters()函数进行简单指数平滑预测
rainseriesforecasts <- HoltWinters(rainseries, beta=FALSE,gamma= FALSE)#其他两个参数我们都设定为FALSE
rainseriesforecasts#返回平滑参数，其实已经确定了公式，当 alpha 越接近 0的时候，临近预测的观测值在预测中的权重就越小

#默认的时候，HoltWinters()默认仅给出在原始时间序列所覆盖时期内的预测,也就是1813-1912的预测值，
rainseriesforecasts$fitted
plot(rainseriesforecasts) #画出预测的值和原有值得比较，发现预测的值更加的平滑
#我们可以计算样本内预测误差的误差平方之和，即原始时间序列覆盖的时期内的预测误差来评估模型的准确度
rainseriesforecasts$SSE#预测误差平方和

#如何预测更多的值呢？使用forecast.HoltWinters(),要先安装包forecast
install.packages("forecast")
library(forecast)
rainseriesforecasts2<-forecast.HoltWinters(rainseriesforecasts,h=8)#预测未来8年的数据
rainseriesforecasts2#同时给出了80%与90%的预测区间
plot.forecast(rainseriesforecasts2)#画出预测值以及预测区间，蓝线是预测值，深色区域是80%预测区间，浅色是90%预测区间

#使用forecast.HoltWinters()返回的样本内预测误差将被存储在一个元素名为“residuals”的列表变量中，如果
#预测模型不可再被优化，连续预测中的预测误差是不相关的。换句话说，如果连续预测中的误差是相关的，很有
#可能是简单指数平滑预测可以被另一种预测技术优化。
#我们使用acf()命令来计算预测误差的相关图，一般计算1-20阶
acf(rainseriesforecasts2$residuals,lag.max=20)#计算1-20阶的预测误差相关图
#可以从样本相关图中看出自相关系数在 3 阶的时候触及了置信界限。为了验证在滞后 1-20 阶（lags 1-20）时
#候的非 0 相关是否显著，我们可以使用 Ljung-Box 检验。
#使用Box.test来计算相关是否显著
Box.test(rainseriesforecasts2$residuals, lag=20,type="Ljung-Box")

# 这里 Ljung-Box 检验统计量为 17.4，并且 P 值是 0.6，所以这是不足以证明样本内预测误差在 1-20阶是非零自相关的。
#为了确定预测模型不可继续优化，我们需要一个好的方法来检验预测误差是正态分布，并且均值为零，方差不变。为了检验预测误差是方差不变的，我们可以画一个样本内预测误差图：
plot.ts(rainseriesforecasts2$residuals)

#为了检验预测误差是均值为零的正太分布，我们可以画出预测误差的直方图，并覆盖上均值为零、标准方差的正
#态分布的曲线图到预测误差上。

#自定义了一个函数
plotForecastErrors <- function(forecasterrors)
{
  # make a red histogram of the forecasterrors:
  mybinsize <- IQR(forecasterrors)/4
  mymin <- min(forecasterrors)*3
  mymax <- max(forecasterrors)*3
  mybins <- seq(mymin, mymax, mybinsize)
  hist(forecasterrors, col="red", freq=FALSE,breaks=mybins)
  # freq=FALSE ensures the area under thehistogram = 1
  mysd <- sd(forecasterrors)
  # generate normally distributed data with mean 0and standard deviation mysd
  mynorm <- rnorm(10000, mean=0, sd=mysd)
  myhist <- hist(mynorm, plot=FALSE,breaks=mybins)
  # plot the normal curve as a blue line on top ofthe histogram of forecast errors:
  points(myhist$mids, myhist$density, type="l",col="blue", lwd=2)
}
plotForecastErrors(rainseriesforecasts2$residuals)

#霍尔特指数平滑法
#如果你的时间序列可以被描述为一个增长或降低趋势的、没有季节性的相加模型，你可以使用霍尔特指数平滑法
#对其进行短期预测。

skirts <-scan("//robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
skirts
skirtsseries<-ts(skirts)
plot.ts(skirtsseries)
skirtsseriesforecasts<-HoltWinters(skirtsseries,gamma=FALSE)
skirtsseriesforecasts
skirtsseriesforecasts$SSE
plot(skirtsseriesforecasts)
skirtsseriesforecasts$fitted
HoltWinters(skirtsseries, gamma=FALSE, l.start=608,b.start=9)
skirtsseriesforecasts2 <-forecast.HoltWinters(skirtsseriesforecasts, h=19)
plot.forecast(skirtsseriesforecasts2)
skirtsseriesforecasts2$fitted
acf(skirtsseriesforecasts2$residuals, lag.max=20)
Box.test(skirtsseriesforecasts2$residuals,lag=20,type="Ljung-Box")
plot(skirtsseriesforecasts2$residuals)
plotForecastErrors(skirtsseriesforecasts2$residuals)


#Holt-Winters 指数平滑法
#如果你有一个增长或降低趋势并存在季节性可被描述成为相加模型的时间序列，你可以使用霍尔特-温特指数平滑
#法对其进行短期预测。
#使用昆士兰海滨度假圣地的纪念品商店从 1987 年 1 月到 1987 年 12 月的每月销售数据
souvenir <-scan("//robjhyndman.com/tsdldata/data/fancy.dat")

souvenirtimeseries<-ts(souvenir,start=c(1987),frequency=12)
plot.ts(souvenirtimeseries)
logsouvenirtimeseries <- log(souvenirtimeseries)
plot.ts(logsouvenirtimeseries)

souvenirtimeseriesforecasts<-HoltWinters(logsouvenirtimeseries)
souvenirtimeseriesforecasts
plot(souvenirtimeseriesforecasts)
souvenirtimeseriesforecasts2<-forecast.HoltWinters(souvenirtimeseriesforecasts,h=48)
plot(souvenirtimeseriesforecasts2)
acf(souvenirtimeseriesforecasts2$residuals,lag.max=20)
Box.test(souvenirtimeseriesforecasts2$residuals,lag=20,type="Ljung-Box")
plot.ts(souvenirtimeseriesforecasts2$residuals)            # make a time plot
plotForecastErrors(souvenirtimeseriesforecasts2$residuals) # make a histogram


#arima模型
#ARIMA 模型为平稳时间序列定义的。因此，如果你从一个非平稳的时间序列开始，首先你就需要做时间序列差分直
#到你得到一个平稳时间序列。如果你必须对时间序列做 d 阶差分才能得到一个平稳序列，那么你就使用
#ARIMA(p,d,q)模型，其中 d 是差分的阶数
#在R中使用diff()函数
#我们使用裙子数据
skirtsseries
plot(skirtsseries)
skirtsseriesdiff<-diff(skirtsseries,differences=1)#一阶差分
plot(skirtsseriesdiff)#一阶差分并不是很平稳
skirtsseriesdiff2<-diff(skirtsseries,differences=2)#再次差分
plot(skirtsseriesdiff2)#二阶差分就好多了

#对于平稳性正式的检验，对于平稳性正式的检验称作“单位根测试”，可以在 fUnitRoots 包中得到，可以在 CRAN 中获得fUnitRoots包而运行

#如果你需要对你的原始时间序列数据做 d 阶差分来获得一个平稳时间序列，那么意味着你可以对你的时间序列使用
#ARIMA(p,d,q,)模型，其中 d 是差分的阶数，裙子数据就是一个arima(p,2,q)的时间序列


#对国王年龄数据进行差分
kingtimeseriesdiff1 <- diff(kingstimeseries,differences=1)
plot(kingtimeseriesdiff1)
#发现一阶差分之后序列就变的平稳了，其实是相当于去除了序列中的趋势部分
#差分之后选择合适的arima模型，其实就是确定p,q的值，这里我们需要用到两个函数acf()，pacf(),建立自相关和偏相关图来看
acf(kingtimeseriesdiff1,lag.max=20)#画自相关图
acf(kingtimeseriesdiff1,lag.max=20,plot=FALSE)#求自相关值

#根据图形确定是1阶截尾，所以是MA(1)

pacf(kingtimeseriesdiff1,lag.max=20)#画偏相关图
pacf(kingtimeseriesdiff1,lag.max=20,plot=FALSE)#求偏相关值
#根据图形确定是3阶截尾，所以是AR(3)
#最终是选择ar还是ma需要考虑aic

#快捷方式：auto.arima()  函数
#我们用国王年龄来进行试验
library(forecast)
kings
auto.arima(kings)#自动确定arima模型的参数，最终是（0，1，1）

#我们用1500-1969 年在北半球火山灰覆盖指数数据（从 Mcleod 和 Hipel得到的原始数据，1994）来做一下arima模型的预测

volcanodust <-scan("//robjhyndman.com/tsdldata/annual/dvi.dat",skip=1)#读数据
volcanodustseries<-ts(volcanodust,start=c(1500))
plot.ts(volcanodustseries)
acf(volcanodustseries,lag.max=20)
pacf(volcanodustseries,lag.max=20)
#确定模型是arma(1,3)

#
auto.arima(volcanodust)#自动确定的模型是arima(1,0,2)

auto.arima(volcanodust,ic="bic")#使用bic来确定参数

#使用arima模型进行预测

kingstimeseriesarima <- arima(kingstimeseries,order=c(0,1,1))  # fit an ARIMA(0,1,1) model
kingstimeseriesarima

library("forecast") # load the "forecast" R library
kingstimeseriesarimaforecasts<-forecast.Arima(kingstimeseriesarima,h=5)#预测5个国王的年龄
kingstimeseriesarimaforecasts
plot.forecast(kingstimeseriesarimaforecasts)

#检验预测误差是否为均值为零，方差为常数的正态分布
acf(kingstimeseriesarimaforecasts$residuals, lag.max=20)
Box.test(kingstimeseriesarimaforecasts$residuals, lag=20,type="Ljung-Box")
plot.ts(kingstimeseriesarimaforecasts$residuals)
plotForecastErrors(kingstimeseriesarimaforecasts$residuals)


#火山灰数据
volcanodustseriesarima<-arima(volcanodust,order=c(1,0,2))
volcanodustseriesarima
library("forecast")
volcanodustseriesforecasts<-forecast.Arima(volcanodustseriesarima,h=31)
plot.forecast(volcanodustseriesforecasts)

#检验预测误差是否为均值为零，方差为常数的正态分布没有证据证明在滞后 1-20 阶（lags1-20）中预测误差是非零自相关的
acf(volcanodustseriesforecasts$residuals, lag.max=20)
Box.test(volcanodustseriesforecasts$residuals, lag=20,type="Ljung-Box")

#检验预测误差是否为均值为零，方差为常数的正态分布
plot.ts(volcanodustseriesforecasts$residuals)
plotForecastErrors(volcanodustseriesforecasts$residuals)