#时间序列

#横截面数据 cross-sectional
#在横截面数据集中，我们是在一个给定的时间点测量变量值

#纵向数据 longitudinal
#随时间的变化饭粗测量变量值
#将研究在给定的一段时间内有规律的记录的观测值，对于这样的观测值，我们可以将其整合成形如Y1,Y2...的时间序列

#对于时序数据基本研究两个问题
#1  对数据的描述，这段时间内发生了什么
#2  预测，接下来会发生什么

#生成时序对象
#time-series object
#即R中包括观测值，起始时间，终止时间以及周期（如月，季度或年）的结构

#范例一
#myseries <- ts(data,start=,end=,frequency=)
#myseries    存储所生成的时序对象
#data        原始的包含观测值的数值型向量
#start       时序的起始时间
#end         时序的终止时间
#frequency   为每个单位时间所包含的观测值数量，1对应年度数据，12对应月度数据，4对应季度数据

sales <- c(18,33,41,7,34,35,24,25,24,21,25,20,
22,31,40,29,25,21,22,54,31,25,26,35) #初始向量类型
tsales <- ts(sales,start=c(2003,1),frequency=12) #强制转换成时序对象
plot(tsales)#画出时间序列的折线图，时间对应横轴，数值对应纵轴
start(tsales) #返回时间序列的开始时间，2003    1
end(tsales)#返回时间序列的结束时间，2004   12
frequency(tsales)#返回时间序列中时间点的个数， 12

#window(),对时序对象取子集
#时序对象，起始时间点，结束时间点
#注意连续持续取值，从起始时间点，到结束时间点
tsales.subset <- window(tsales,start=c(2003,5),end=c(2004,6))

#第一步 画图

#时序数据集中通常有很显著的随机或误差成分。为了辨明数据中的规律，我们总是希望
#能够撇开这些波动，画出一条平滑曲线
#最简单的方法是简单移动平均

install.packages("forecast")
install.packages("xts")
library(zoo)
library(xts)

library(forecast)

opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
ylim <- c(min(Nile),max(Nile))

#ma(),拟合一个简单的移动平均模型
#并没有什么特别的科学理论来指导k的选取，只是需要先尝试多个不同的k,再决定一个最好的k
#从图像来看，随着k的增大，图像会变的越来越平滑
#需要找到最能画出规律的k,避免过平滑或者欠平滑
plot(ma(Nile,3),main="Simple Moving Averages (k=3)",ylim=ylim)
plot(ma(Nile,7),main="Simple Moving Averages (k=7)",ylim=ylim)
plot(ma(Nile,15),main="Simple Moving Averages (k=15)",ylim=ylim)
par(opar)

#对于间隔大于1的时序数据（即存在季节性因子），我们需要了解的就不仅仅是总体趋势了
#此时我们需要通过季节性分解来帮忙我们探究季节性波动以及总体趋势

#季节性分解
#存在季节性因素的时间序列数据(如月度数据，季度数据等)可以被分解为：
#趋势因子，trend component,能捕捉到长期变化
#季节性因子，seasonal component能捕捉到一年内的周期性变化
#随机(误差)因子,irregular/error component,则能捕捉到那些不能被趋势或季节效应解释的变化

#利用模型来分解数据

#相加模型
#各种因子之和应等于对应的时序值
#y=Trend+Seasonal+Irregular

#相乘模型
#各个因子之积等于对应的时序值
#y=Trend*Seasonal*Irregular

#比较相加模型和相乘模型
#假设有一个时序，记录了10年来摩托车的月销售量
#可加模型中，11月，12月(圣诞节)的销量一般会增加500，而1月的销量则会减少200
#此时季节性的波动量和当时的销量无关
#相乘模型中,11月，12月的销量则会增加20%，1月的销量减少10%，
#即季节性的波动量和当时的销量是成比例的。
#很多时候，相乘模型比相加模型更现实一些

#范例一
plot(AirPassengers)  #大值数据直接画图
lAirPassengers <- log(AirPassengers)#将大值进行对数转换，转换成小一号的值5-7之间
plot(lAirPassengers,ylab="log(AirPassengers")#转换后的对数数据直接画图

#用LOESS光滑将时序分解成季节项，趋势项和随机项
#stl只能处理相加模型，但这也不算一个多严重的限制，因为相乘模型总可以通过对数变换
#转换成相加模型，用经过对数变换的序列拟合出的相加模型也总可以再转换回原始尺寸

#stl(ts,s.windaow=,t.window=)
#ts,        将分解的时序
#s.window   控制季节效应变化的速度,periodic,可使得季节效应在各年间都一样
#t.window   控制趋势项变化的速度
#较小的值意味着更快的变化速度

#原始对数时序行列中的每一个值，变成行标识，每一个值都分解成3部分
fit <- stl(lAirPassengers,s.window="period")
class(fit) #返回"stl"
plot(fit)
#图形解读
#data       图形第一块，对应时序图
#seasonal   图形第二块，对应季节效应图
#trend      图形第三块，对应趋势图
#remainder  图形第四块，对应随机波动图
#序列图的趋势是单调增长，季节效应表明夏季乘客数量更多(可能因为假期)
#每个图的y轴尺度不同，因此我们通过图中古侧的灰色长条来指示量级，即每个长条代表的量级一样

#具体参看每个观测值各分解项的值
fit$time.series  #取过Log

#指数函数，返回到取对数前的值,方便参数解释
exp(fit$time.series)
#参数解释
#         seasonal    trend     remainder
#Jul 1957 1.2415994  366.2637   1.0225336
#seasonal=1.2415994
#Jul 对应7月
#说明7月的乘客数增加了24%（0.24=1.24-1）

#         seasonal    trend     remainder
#Nov 1959 0.8077298 442.3741 1.0131012
#seasonal=0.8077298
#Nov 对应11月
#说明11月的乘客数减少了20%（0.2=0.8077298-1)

#monthplot 对季节分解进行可视化
par(mfrow=c(2,1))
library(forecast)
#月度图
#表示的是每个月份组成的子序列，
#J~D,表示1月=12月，J上的一根折线，表示1月份的子序列，包含所有年份中的1月份
#连接所有1月的点，连接2月所有的点，依次类推，以及每个子序列的平均值
#从这幅图看，每个月的增长趋势几乎是一致的
monthplot(AirPassengers,xlab="",ylab="")
#季节图
#这幅图以年份为子序列
#从这幅图看，每一年当中7月最高，每一年趋势都基本相同，与月度图类似
seasonplot(AirPassengers,year.labels="TRUE",main="")

#指数预测模型
#指数模型是用来预测时序未来值的最常用模型

#范例一 单指数平滑
#单指数模型 simple/single exponential model
#拟合的是只有常数水平项和时间点i处随即向的时间序列，这时认为时间序列不存在趋势项和季节项
#根据现有的时序值的加权平均对未来值做短期预测，其中权数选择的宗旨是使得距离现在越远的观测值
#对平均数的影响越小

library(forecast)
#ets(ts,model="ZZZ)
#ts,要分析的时序
#model，限定模型的字母
#位置说明：
#第一个字母，代表误差项
#第二个字母，代表趋势项
#第三个字母，代表季节项
#相加模型 A
#相乘模型 M
#无       N
#自动选择 Z

#尝试对时间序列对象nhtemp进行季节分组
stl(nhtemp,s.window = "period")
#报错：Error in stl(nhtemp, s.window = "period") : 序列没有周期，或其周期小于二
#查看时间序列对象 Time Series
#Time Series:
#Start = 1912
#End = 1971
#Frequency = 1
#本身单周期，不需要进行季节分组，即仅1列数据，无法进行季节性分解

fit <- ets(nhtemp,model="ANN") #对nhtemp时序对象拟合单指数模型
fit
#结果解读：
#ANN，误差项相加模型，趋势项模型为空，季节项模型为空
#alpha,平滑指数，控制水平向的指数型下降
#为拟合某种标准，实际值一般由计算机选择，常见的拟合标准是真实值和预测值之间的残差平方和
#alpha = 0.1819 ，说明预测同时考虑了离现在较近和较远的观测值，这样的值可以最优化模型在给定数据集上的拟合效果

#forecast(fit,k)，用于预测时序未来的k步
forecast(fit,1)
#结果解读:
#Point Forecast=51.87031,返回执行一步向前预测后的时序值
#Lo 80    Hi 80
#50.40226 53.33835  80%的置信区间[50.4,53.3]
#Lo 95    Hi 95
#49.62512 54.11549  95%的置信区间[49.7,54.1]

#画原始时序图
plot(nhtemp)
#画预测时序图
#结果解读：
#相当于在原始时序图末尾，增加一个预测，纵坐标为时序值，标注了预测值80%，以及95%的置信区间
plot(forecast(fit,1),xlab="Year",
ylab=expression(paste("Temperature (",degree*F,")")),
main="New Haven Annual Mean Temperature")

#返回时序模型的拟合优度度量
accuracy(fit)
#结果解释
#ME       平均误差mean(e)，正向和负向的误差会相互抵消，用处不大
#RMSE     平均残差平方和的平方根 ，最有名，最常用

#MAE      平均绝对误差
#平均绝对误差给出了误差在真实值中的占比，它没有单位，因此可以用于比较不同时序间的预测准确性，
#但它同时假定测量尺度中存在一个真实为零的值，比如每天的游客数量，但华氏温度中并没有一个真实为零
#的值，因此这里不能用这个统计量。

#MPE      平均百分比误差，正向和负向的误差会相互抵消，用处不大
#MAPE     平均绝对百分误差

#MASE     平均绝对标准化误差
#最新的一种准确度测量
#通常用于比较不同尺度的时序间的预测准确性

#双指数平滑
#双指数模型 double exponential model 也叫Holt指数平滑，Holt exponential smoothing
#Holt指数平滑，拟合了有水平项和趋势项的时序

library(forecast)

#查看时序数据源基本信息
class(AirPassengers) #返回ts
frequency(AirPassengers) #返回时间序列中时间点的个数，即列个数12，满足季节性分解的条件

#需要对原始数据取对数，使得其满足可加模型
fit <- ets(log(AirPassengers),model="AAA")
class(fit)  #返回ets对象
fit
#结果解读：
#ETS(A,A,A) 对应model="AAA"

#alpha = 0.6975，平滑参数alpha，控制水平项指数型下降，水平项
#beta  = 0.0031，平滑参数beta,  控制斜率指数型下降，趋势项
#gamma = 1e-04 ，平滑参数gamma, 控制季节项的指数下降，季节项
#三个参数的有效范围都是[0,1],参数值越大意味着观测值的权重越大
#趋势项的参数小意味着近期观测值的斜率不需要更新

accuracy(fit)
#结果解读：
#展示了时序预测中的几个最主流的几个准确性度量

pred <- forecast(fit,5)
pred
#结果解读:
#5步向前预测的结果是对应5行值，一步向前就对应一行值

#画原时序图
plot(log(AirPassengers))
#画预测时序图
plot(pred,main="Forecast for Air Travle",
ylab="Log(AirPassengers)",xlab="Time")
#结果解读：
#将后续预测的5个月份的时序值，以及对应的80%,95%的置信区间续貂到原图尾部

pred$mean <- exp(pred$mean) #时序值，从对数逆运算，返回到指数运算
pred$lower <- exp(pred$lower) #80%，95%置信区间下端，从对数逆运算，返回到指数运算
pred$upper <- exp(pred$upper)#80%，95%置信区间上端， 从对数逆运算，返回到指数运算
p <- cbind(pred$mean,pred$lower,pred$upper)#基于向量，创建数据框

#dimnames(p) 返回数据框行名，列名
#class(dimnames(p)) #返回列表
#[[1]]
#NULL

#[[2]]
#[1] "mean"  "Lo 80" "Lo 95" "Hi 80" "Hi 95"

dimnames(p)[[2]] <- c("mean","Lo 80","Lo 95","Hi 80","Hi 95") #给列名赋值
p#改造pred为指数运算后的pred
#结果解读
#Mar 1961 511.1312 477.0088 459.8775 547.6945 568.0971
#三月份，将有512个乘客，95%的置信区间[459.8775,568.0971]

#ets函数和自动预测

#范例一 使用ets 进行自动指数预测
library(forecast)

#查看数据源
class(JohnsonJohnson)  #返回ts
frequency(JohnsonJohnson) #返回时间序列中时间点的个数，即对应列数=4，满足季节性分解

#ets
#还可以用来拟合有可乘项的指数模型
#ets(log(AirPassengers),model="AAA"),对数，可加模型
#也支持ets(AirPassengers,model="MAM")，原始，可乘模型
#当前，仍假定趋势项可加，季节项，误差项可乘

#可以用来拟合抑制项
#时序预测一般假定序列的长期趋势是一直向上(如房价市场)，
#而一个抑制项则使得趋势项在一段时间内靠近一条水平渐近线
#在许多问题中，一个有抑制项的模型往往更符合实际情况

#不指定model时，自动选取对原始数据拟合优度最高的模型
fit <- ets(JohnsonJohnson)
fit
#结果解读：
#ETS(M,A,A) 对应模型MAA,
#水平项(或者称为误差项)相乘，趋势项相加，季节项相加

#画原图
plot(JohnsonJohnson)
#画预测图
#forecast(fit,k)，k缺省为8
plot(forecast(fit),main="Johnson & Johnson Forecasts",
ylab="Quarterly Earnings (Dollars)",
xlab="Time",flty=2)

#ARIMA预测模型
#在ARIMA预测模型中，预测值表示为最近的真实值和最近的预测误差组成的线性函数

#建立ARIMA模型的步骤包括
#第一步：确保时序是平稳的
#第二步：找到一个(或几个)合理的模型(即选定可能的p值和q值)
#第三步：拟合模型
#第四步：从统计假设和预测准确性等角度评估模型
#第五步：预测

#范例一 序列的变换以及稳定性评估 第一步
library(forecast)
library(tseries)

#查看数据源
#class(Nile)  #返回ts类型
#frequency(Nile) #返回1，列数=1，不能进行季节性分解

#画原图
#平稳性一般可以通过时序图直观判断。
#如果方差不是常数，则需要对数据做变换。
#如果数据存在趋势项，则需要对其进行差分
plot(Nile)
#结果解读：
#各观测年间的方差似乎是稳定的，因此我们无需对数据做变换
#但数据中可能存在某种下降趋势

#找到最优的d值
#由于一般假定平稳性时序有常数均值，这样的序列肯定不含有趋势项。
#非平稳的时序可以通过差分类转换成平稳性序列
#具体来说，差分就是将时序中每一个观测值Yt都替换成Y(t-1)-Yt
#注意对序列的一次差分可以移除序列中的线性趋势
#二次差分移除二次项趋势，三次差分移除三次项趋势
#在实际操作中，对序列进行两次以上的差分通常都是不必要的
ndiffs(Nile)  #返回1，即最优d值=1

#diff(ts,differences=d)
#其中d即对序列ts的差分次数，默认值d=1
dNile <-diff(Nile)
#结果解读：
#class(dNile)  #返回ts

#画经过差分以后的图
plot(dNile)
#结果解读：
#差分后原始数据中下降的趋势被移除掉了

#ADF(Augmented Dickey-Fuller)统计检验来验证平稳性假设
adf.test(dNile)
#p-value== 0.01,落入拒绝域，接受原假设，序列是平稳的

#范例二  选择模型   第二步
#选择ARIMA模型的方法,在滞后阶数后面的对应
# 模型                ACF                    PACF
#ARIMA(p,d,0),    逐渐减少到零          在p阶后减少到零
#ARIMA(0,d,q),    在q阶后减少到零       逐渐减少到零
#ARIMA(p,d,q),    逐渐减少到零          逐渐减少到零

par(mfrow=c(2,1))

#一阶滞后(Lag 1),代表时序向左移动一位
Nile #Time Series:Start = 1871 End = 1970
lag(Nile,0) #Time Series:Start = 1871 End = 1970,等于原始值
lag(Nile,1) #Time Series:Start = 1870 End = 1969,所有值不动，时间轴向左移动一位，减1
#AC1就是一阶滞后序列和0阶滞后序列间的相关性
#(AC1,AC2....,ACk)构成的图即为自相关函数图
#自相关函数图 AutoCorrelation Function plot
Acf(dNile) #生成Acf图

#偏自相关
#即当序列Yt和Y(t-1)之间的所有值带来的效应都被移除后，两个序列间的相关性
Pacf(dNile)#生成Pacf图

#结果解释:
#ACF,PACF图形的横坐标对应滞后阶数
#参考线为0
#ACF图中，滞后阶数为1时，有一个比较明显的自相关，当lag增加，自相关逐渐减少至0
#PACF图中，滞后阶数为1时，而当其继续增加时，偏相关逐渐减少至零
#因此参考选择ARIMA模型的方法第二行：
#其中ndiffs(Nile)  #返回1，即最优d值=1
#滞后阶数为1，表明q=1
#因此，建议考虑ARIMA(0,1,1)模型

#范例三  拟合ARIMA模型  第三步
library(forecast)
fit <- arima(Nile,order=c(0,1,1))
fit
#结果解读;
#模型为ARIMA(0,1,1)
#移动平均项的系数 -0.7329
#模型的AIC=1269.09,
#如果我们还有其他备选模型，则可以通过比较AIC值来得到最合理的模型
#AIC值越小的模型越好

accuracy(fit)

#范例四   模型评价  第四步
#一般来说，一个模型如果合适，那模型的残差应该满足均值为0的正态分布
#并且对于任意的滞后阶数，残差自相关系数都应该为零
#换句话说，模型的残差应该满足独立正态分布（即残差间没有关联）

#如果满足正态分布，则数据中的点会落在图中的线上
qqnorm(fit$residuals)
qqline(fit$residuals)

#Box.test 函数可以检验残差的自相关系数是否都为零
Box.test(fit$residuals,type="Ljung-Box")
#解读结果：
#原假设：残差的自相关系数=0
#p-value = 0.2416，未落入拒绝域，接受原假设
#即认为残差的自相关系数为0
#如果模型残差不满足正态性假设或零自相关系数假设，则需要调整模型，
#增加参数或改变差分次数

#范例五 预测  第五步
#原始Nile #Time Series:Start = 1871 End = 1970
#利用fit,对接下来的三年进行预测，1971,1972,1973
forecast(fit,3)
plot(forecast(fit,3),xlab="Year",ylab="Annual Flow")

par(mfrow=c(1,1))

#ARIMA的自动预测

#范例
library(forecast)

#观察数据源
class(sunspots)      #返回ts类型
frequency(sunspots)  #返回12，可以做季节性分解
start(sunspots)      #返回1749 1，时序对象的起始时间
end(sunspots)        #返回1983 12，时序对象的截止时间

fit <- auto.arima(sunspots)
fit
#结果解读：
#最后选中模型ARIMA(2,1,2)，与其他模型相比，在这种情况下得到的模型的AIC值最小
#

accuracy(fit)
#结果解读：
#MPE=Nan,MAPE=Inf
#由于观测值中存在值为0的观测，所有MPE,MAPE这两个度量值失效

#end(sunspots)        #返回1983 12，时序对象的截止时间
#在原始截止时间的基础上，预测后3个数据
#Jan 1984 ，Feb 1984，Mar 1984
forecast(fit,3)

#时间序列#横截面数据 cross-sectional#在横截面数据集中，我们是在一个给定的时间点测量变量值#纵向数据 longitudinal#随时间的变化饭粗测量变量值#将研究在给定的一段时间内有规律的记录的观测值，对于这样的观测值，我们可以将其整合成形如Y1,Y2...的时间序列#对于时序数据基本研究两个问题#1  对数据的描述，这段时间内发生了什么#2  预测，接下来... 　　模式发生时间的一个有序 序列 ，se=<t1,t2,t3,t4,t5> 2、 周期 点 序列 <0,5,10,15,20,27,30,35,40>是一个 周期 为5、时间容忍度为2的 周期 点 序列 。 3、部分 周期 模式 　　<0,5,10,15,27,30,35,40>可以分为<0,5,10,15>、<15,27>和<27... dat <- read_excel("rawdatacompet.xlsx",sheet=1,na="NA") series1 <- dat[,1] series1 <- series1[1:60,] #将数据转化为时间 序列 series1 <- ts(series1) #将数据转化为双精度， 没有 貌似 时间 序列 提供了预测未来价值的机会。基于以前的价值观，可以使用时间 序列 来预测经济，天气和能力规划的趋势，其中仅举几例。时间 序列 数据的具体属性意味着通常需要专门的统计方法。