网站搭建响应式,阜阳建设网站公司,做外贸女装有哪些网站有哪些,网站制作长春1 数据集说明
AirPassengers 1949~1960年每月乘坐飞机的乘客数
JohnsonJohnson JohnsonJohnson每股季度收入
nhtemp 康涅狄格州纽黑文地区从1912年至1971年每年的平均气温
Nile 尼罗河的流量
sunspots 1749年~1983年月平均太阳黑子数
2 相关包
xts、forecast、tser…1 数据集说明
AirPassengers 1949~1960年每月乘坐飞机的乘客数
JohnsonJohnson JohnsonJohnson每股季度收入
nhtemp 康涅狄格州纽黑文地区从1912年至1971年每年的平均气温
Nile 尼罗河的流量
sunspots 1749年~1983年月平均太阳黑子数
2 相关包
xts、forecast、tseries、directlabels
3 在R中生成时序对象
首先将分析对象转成时间序列对象time-series object。
xts包提供的xts类支持时间间隔规律的和时间间隔不规律的时间序列还拥有很多用于操作时序数据的函数。
另外基础R附带的ts用于储存时间间隔规律的耽搁时间序列mts’用于储存时间间隔均匀的多个时间序列zoo包提供的类可以储存时间间隔不规律的时间序列xts包提供zoo类的超集包含更多函数tsbox包提供的函数可以将数据框转化为任意一种时间序列格式也可以将一种时间序列格式转化为另一种时间序列格式。
创建xts时间序列
library(xts)myseries-xts(data,index)
# data是观测值的数值型向量index表示观测值观测时间的日期向量。
3.1 生成时间序列对象
# 从2018年1月以来的两年的月度销售数据library(xts)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)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)date-seq(fromas.Date(2018/1/1),toas.Date(2019/12/1),bymonth)sales.xts-xts(sales,date)sales.xts
# 结果 xts格式的时间序列对象可使用方括号[]来设定子集
# 返回2018年以来的所有数据sales.xts[2018]
# 结果 # 返回2018年3月到2019年5月的所有数据sales.xts[2018-3/2019-5] 函数apply用于对时间序列对象的每个不同时间段执行一个函数在将时间序列聚集合成更大的时间段时很实用。
newseries-apply.period(x,FUN,…)
# period可以是daily、weekly、monthly、quarterly或yearlyx是一个xts时间序列对象FUN是要应用的函数…是传递给FUN的参数。
# 返回8个季度销售总量的时间序列quarterlies-apply.quarterly(sales.xts,sum)quarterlies
# 结果 # sum可以替换为mean、median、min、max或其他任何返回单个值的函数
3.2 绘制时间序列图
forecast包中的autoplot()可将时序数据绘制为ggplot2图形
library(ggplot2)library(forecast)autoplot(sales.xts)
# 结果 图 3-1 销量数据的时序图autoplot()提供的默认格式
autoplot(sales.xts)geom_line(colorblue)scale_x_date(date_breaks1 months,date_labels%b %y)labs(x,ySales,titleCustomized Time Series Plot)theme_bw()theme(axis.text.xelement_text(angle90,vjust0.5,hjust1),panel.grid.minor.xelement_blank())
# 选项date_breaks指定刻度间的距离值可以是“1 day”“2 weeks”“5 years”等“%b%y”指定了月份3个字母和年份2位数以及两者之间的空格最后选择了黑白主题x轴标签旋转了90度取消了垂直的小网格线。
# 结果 图 3-2 销量数据时序图定义了颜色、更美观的标签和更清晰的主题元素
4 时序的平滑化和季节项分解
对时序数据建立复杂模型之前需要对其进行描述和可视化。我们对时序进行平滑化以探究其总体趋势并对其进行分解以观察时序中是否存在季节项。
4.1 通过简单移动平均进行平滑处理Nile数据集
处理时序数据的第一步是画图。使用数据集Nile是埃及阿斯旺市在1871年至1970年间所记录的尼罗河的年度流量。由图数据总体呈下降趋势但不同年份的变动非常大。图 4-1
时序数据集中通常有很显著的随机或误差成分。为了辨明数据中的规律我们希望能够撇开这些波动画出一条平滑曲线。画出平滑曲线的最简单方法是简单移动平均。比如每个数据点都可用这一点和其前后两个点的平均值来表示这就是居中移动平均。
R中可以做简单移动平均的函数TTR包中的SMA()函数zoo包中的rollmean()函数forecast包中的ma()函数。
# 时序数据的原始数据图、平滑后的图library(ggplot2)library(forecast)theme_set(theme_bw())ylim-c(min(Nile),max(Nile))autoplot(Nile)ggtitle(“Raw time series”)scale_y_continuous(limitsylim)autoplot(ma(Nile,3))ggtitle(“Simple Moving Average (k3)”)scale_y_continuous(limitsylim)autoplot(ma(Nile,7))ggtitle(“simple Moving Average (k7)”)scale_ycontinuous(limitsylim)autoplot(ma(Nile,15))ggtitle(“simple Moving Average (k15)”)scale_ycontinuous(limitsylim)
# 结果 图 4-1 1871~1970年阿斯旺水站观测到的尼罗河的年流量 图 4-2 平滑水平k3 图 4-3 平滑水平k7 图 4-4 平滑水平k15
# 随着k的增大图像变得越来越平滑。
# 尼罗河的·流量从1892年到1900年有明显下降其他的变动不是太好解读比如1941年到1961年流量似乎略有上升但这也可能只是一个随机波动。
我们需要找到最能画出数据中规律的k避免过平滑或者欠平滑。没有什么特别的科学理论来指导k的选取只需要多次尝试不同的k再决定一个最好的k。
对于间隔大于1的时序数据存在季节想需要了解的不仅是总体趋势。此时需要通过季节项分解帮助我们探究季节性波动以及总体趋势。
4.2 季节项分解
存在季节项的时序数据如月度数据、季度数据等可以被分解为趋势项、季节项和随机项。趋势项trend component能捕捉到长期变化季节项seasonal component能捕捉到一年内的周期性变化随机误差项irregular/error component能捕捉到那些不能被趋势项或季节项解释的变化。
可以通过相加模型或相乘模型来分解数据。在相加模型中各项之和等于对应的时序值时刻t的观测值等于这一时刻的趋势项、季节项以及随机项之和而在相乘模型中时刻t的观测值等于趋势项、季节项和随机项相乘。
下面举例说明相加模型与相乘模型的区别假设有一个时序记录了10年来摩托车的月销量。在具有季节项的相加模型中11月和12月的销量一般会增加500圣诞节销售高峰而1月一般是销售淡季的销量则会减少200.此时季节性的波动量和当时的销量无关。在具有季节项的相乘模型中11月和12月的销量会增加20%1月的销量会减少10%即季节项的波动量和当时的销量是成比例的这与相加模型不一样。这也使得在很多时候相乘模型比相加模型更切合实际。
将时序分解为趋势项、季节项和随机项的常用方法是用loess平滑做季节项分解通过stl()函数实现
stl(ts, s.window, t.window)
# ts是将要分解的时序s.window控制季节项变化的速度t.window控制趋势项变化的速度。设置s.windows”periodic”可使得季节项宰割年检都一样。参数ts和s.windows必须提供。
# stl()函数只能处理相加模型但是相乘模型可以通过对数变换转换成相加模型。
例 使用数据集AirPassengers
该时序数据集描述了1949年~1960年每个月国际航班的乘客数单位千。
library(ggplot2)library(forecast)autoplot(AirPassengers)
# 结果 图 4-5 AirPassengers时序图
# 序列的波动随着整体水平的增长而增长相乘模型更适合这个序列。
# 对数变换AirPassengers1-log(AirPassengers)autoplot(AirPassengers1,ylablog(AirPassengers))
# 结果 # 分解时间序列## 季节项分解将季节项限定为每年都一样fit-stl(AirPassengers1,s.window”period”)autoplot(fit)
# 结果 图 时序图、季节项分解图、趋势项分解图、随机项分解图
# 时序的趋势为单调增长季节项表明可能为假期夏季乘客更多。
# 每个图的y轴尺度不同通过图中右侧的灰色长条来指示量级即每个长条代表的量级一样
# 返回对数变换后的时序fit$time.series
# 结果部分 # stl()函数返回的对象中有一项是time.series包括每个观测值中的各分解项——趋势项、季节项以及随机项的值。
# 将结果转化为原始尺度exp(fit$time$series)
# 结果部分 # 观察季节项7月的乘客数增长了24%乘子为1.2411月的乘客数减少了20%乘子为0.8
例 通过forecast包提供的其它工具对季节项分解进行可视化月度图和季节图
library(forecast)library(ggplot2)library(directlabels)# 月度图ggmonthplot(AirPassengers)labs(title”Month plot: AirPassengers”,x””,y”Passengers (thousands)”)
# 结果 图 AirPassengers时序的月度图。月度图显示了按月划分的子序列从1949年到1960年每年的所有1月的点连接起来所有2月的点连接起来以此类推以及每个字序列的平均值。每个月的增长趋势几乎一致大多数乘客在7月和8月乘坐飞机。
# 季节图p-ggseasonplot(AirPassengers)geom_point()labs(title”Seasonal plot: AirPassengers”,x””,y”Passengers (thousands)”)direct.label(p)
# 结果 图 季节图显示每年的子序列。我们可以从图中看到类似的规律包括每年乘客的增长趋势和相同的季节性规模。默认情况下ggplot2包将为年份变量创建图例。使用directlabels包可以将年份标签直接放置在图中时序的每条线旁边。
5 指数预测模型
指数模型是用来预测时序未来值的最常用模型。这类模型相对比较简单实践证明它们的短期预测能力良好。不同指数模型建模时选用的因子可能不同。比如单指数模型simple/single exponential model拟合的是只有水平项和时间点i处随机项的时间序列这时认为时间序列不存在趋势项和季节项双指数模型double exponential model也叫Holt指数平滑Holt exponential smoothing拟合的是有水平项和趋势项的时序三指数模型triple exponential model也叫Holt-Winters指数平滑Holt-Winters exponential smoothing拟合的是有水平项、趋势项以及季节项的时序。
Forecast包中的ets()函数可以拟合指数模型
ets(ts,model”zzz”)
# ts是要分析的时序限定模型的字母有3个第1个字母代表随机项第2个字母代表趋势项第3个字母代表季节项。可选的字母包括A相加模型、M相乘模型、N无、z自动选择。
用于拟合单指数、双指数和三指数预测模型的函数 类型 参数 函数 单指数 水平项 ets(ts,model”ANN”) ses(ts) 双指数 水平项、趋势项 Ets(ts,model”AAN”) holt(ts) 三指数 水平项、趋势项、季节项 ets(ts,model”AAA”) hw(ts)
函数ses()、holt()、hw()都是函数ets()的便捷包装函数中有事先默认设定的参数值。
5.1 单指数平滑
单指数平滑根据现有的时序值的加权平均对未来值做短期预测权数选择的宗旨是使得距离现在越远的观测值对平均数的影响越小。
一步向前预测可看做当前值和全部历史值的加权平均。
例 使用数据集nhtemp
nhtemp时序中有康涅狄格州纽黑文市从1912年至1971年每年的平均气温。
# 拟合模型library(forecast)fit-ets(nhtemp,model”ANN”)fit
# A表示可加误差NN表示时序中不存在趋势项和季节项。
# 结果 # α值比较小α0.18说明预测时同时考虑了离现在较近和较远的观测值这样的α值可以最优化模型在给定数据集上的拟合效果。
# 一步向前预测forecast(fit,1)
# 结果 # 一步向前预测的结果是51.87031℉其95%的置信区间为49.62512℉到54.11549℉。
# 输出准确度度量accuracy(fit)
# 结果 # 预测准确度度量 一般来说平均误差和平均百分比误差用处不大因为正向和负向的误差会抵消掉
RMSE给出了平均误差平方和的平方根即1.126℉平均绝对百分比误差给出了误差在真实值中的占比它没有单位可以用于比较不同时序间的预测准确度但它同时假定测量尺度中存在一个真实为零的点比如每天的乘客数但温度中并没有一个真实为零的点因此这里不能使用这个统计量平均绝对标准化误差是最新的一种准确度评估指标通常用于比较不同尺度的时序间的预测准确度。这几种预测准确度评估指标中并不存在某种最优评估指标不过相对来说RMSE最有名、最常用。
5.2 Holt指数平滑和Holt-Winters指数平滑
Holt指数平滑可以对有水平项和趋势项斜率的时序进行拟合。
Ets函数中的平滑参数alpha控制水平项的指数型下降beta控制趋势项的指数型下降。两个参数的有效范围都是[0,1]参数取值越大意味着越近的观测值的权重越大。
Holt-Winters指数平滑可用来拟合含有水平项、趋势项以及季节项的时间序列。此时模型可表示为Ytlevelslope*tstirregulart
st代表时刻t的季节项ets()函数的参数gamma控制季节项的指数下降取值范围是[0,1]取值越大意味着越近的观测值的季节项权重越大。
例 用指数模型预测未来的乘客数Holt-Winters指数平滑
# 有水平项、趋势项以及季节项的指数平滑模型## 平滑参数llibrary(forecast)fit-ets(log(AirPassengers),model”AAA”)fit
# 结果 # 给出了3个平滑参数即水平项0.6795、趋势项0.0031、季节项1e-04。趋势项的值很小意味着近期观测值的趋势不需要更新
accuracy(fit)
# 结果 ## 未来值预测pred-forecast(fit,5)pred
# 结果 autoplot(pred)labs(title”Forecast for Air Travel”,y”Log(AirPassengers)”,x”Time”)
# 结果 基于Holt-Winters指数平滑模型的5年航线乘客数预测对数变换后单位千
## 用原始尺度预测pred$mean-exp(pred$mean)pred$lower-exp(pred$lower)pred$upper-exp(pred$upper)p-cbind(pred$mean,pred$lower,pred$upper)dimnames(p)[[2]]-c(mean,Lo 80,Hi 80,Lo95,Hi 95)p
# 结果 5.3 ets()函数和自动预测
ets()函数还可以用来拟合有可乘项的指数模型加入抑制项dampening component以及进行自动检测。
通过ets(AirPassengers,model”MAM”)函数对原始数据拟合可乘模型。此时仍假定趋势项可加但季节项和误差项可乘。当采用可乘模型时准确度统计量和预测值都基于原始尺度以千为单位的乘客数这也是它的一个明显优势。
ets()函数也可以用来拟合抑制项。时序预测一般假定序列的长期趋势是一直向上的如房地产市场而一个抑制项则使得趋势项在一段时间内靠近一条水平渐近线。在很多问题中一个有抑制项的模型往往更符合实际预测。
最后可以通过ets()函数最懂选取对原始数据拟合优度最高的模型。
例 自动选取最优拟合模型使用ets()函数进行自动指数预测使用数据集JohnsonJohnson
library(forecast)fit-ets(JohnsonJohnson)fit
# 结果 # 没有指定模型软件自动搜索了一系列模型并在其中找到最小化拟合标准默认为对数似然的模型。所选模型同时有可乘趋势项、季节项和误差项随机项。
# 下8个季度的预测autoplot(forecast(fit))labs(x”Time”,y”Quarterly Earnings (Dollars)”,title”Johnson and Johnson Forecasts”)
# 结果 带趋势项和季节项的可乘指数平滑模型其中预测值由虚线表示80%和95%置信区间分别由淡蓝色和深蓝色表示
6 ARIMA预测模型
在ARIMA预测模型中预测值表示为由最近的真实值和最近的预测误差残差组成的线性函数。
6.1 概念介绍
时序的滞后阶数向后追溯的观测值的数量
表 Nile时序的不同滞后阶数 0阶滞后项Lag 0代表没有移位的时序一阶滞后Lag 1代表时序向左移动一位二阶滞后Lag 2代表时序向左移动两位以此类推。
Lag(ts,k)
# 把时序变成k阶滞后ts指代目标序列k为滞后项阶数。
自相关度量时序中各个观测值之间的相关性。相关性构成的图即自相关函数图autocorrelation function plotACF图。ACF图可用于为ARIMA模型选择合适的参数并评估最终模型的拟合效果。
foecast包中的Acf()函数可以生成ACF图
acf(ts)
# ts是原始时序
偏自相关图PACF使用forecast包中的Pacf()函数绘制
Pacf(ts)
PACF图也可以用来找到ARIMA模型中最适宜的参数
ARIMA模型主要用于拟合具有平稳性或可以被转换为平稳序列的时间序列。在一个平稳的时序中序列的统计性质并不会随着时间的推移而改变另外对任意滞后阶数k序列的自相关性不改变。
一般来说拟合ARIMA模型前需要变换序列的值以保证方差为常数比如对数变换、Box-Cox变换等。
一般假定平稳性时序有常数均值这样的序列中肯定不含有趋势项。非平稳的时序可以通过差分来转换为平稳性序列。对序列的一次差分可以移除序列中的线性趋势二次差分移除二次项趋势三次差分移除三次项趋势。在实际操作中对序列进行两次以上的差分通常都是不必要的。
通过函数diff()对序列进行差分
diff(ts,differencesd)
# d是对序列ts的差分次数默认值为d1
# forecast包中的ndiffs()函数可以帮我们找到最优的d值ndiffs(ts)
平稳性一般可以通过时序图直观判断。如果方差不是常数需要对数据做变换如果数据中存在趋势项需要对其进行差分。也可以通过ADFAugmented Dickey-Fuller统计检验来验证平稳性假定。tseries包中的adf.test()函数可以用来做ADF检验adf.test(ts)。如果结果显著则认为序列满足平稳性。
总之通过ACF图和PACF图来为ARIMA模型选定参数。平稳性是ARIMA模型中的一个重要假设可以通过数据变换和差分使得序列满足平稳性嘉定。
下面我们就可以拟合出有自回归autoregressiveAR项、移动平均moving averagesMA项或者两者都有ARMA的模型了最后检验有ARMA项的ARIMA模型并对其进行差分以保证平稳性。
6.2 ARMA和ARIMA模型
ARIMA(p,d,q)模型时序被差分了d次且序列中的每个观测值都是用过去的p个观测值和q个残差的线性组合表示的。预测是“无误差的”或完整的。
建立ARIMA模型的步骤
1确保时序是平稳的
2找到一个或几个合理的模型即选定可能的p值和q值
3拟合模型
4从统计假设和预测准确度等角度评估模型的拟合效果
5预测’
例 对Nile时序拟合ARIMA模型
1验证序列的平稳性
画出序列的折线图并判别其平稳性
library(forecast)library(tseries)autoplot(Nile)
# 结果 1871年~1970年在阿斯旺地区测量的尼罗河的年流量时序图
ndiffs(Nile)
# 结果 dNile-diff(Nile)autoplot(dNile)
# 结果 被差分一次后的时序图差分后原始图中下降的趋势被移除了
adf.test(dNile)
# 结果 # 对差分后的时序做ADF检验结果显示时序此时是平稳的可以继续进行下一步
2选择模型
通过ACF图和PACF图来选择备选模型
autoplot(Acf(dNile))autoplot(Pacf(dNile))
# 结果 # 需要为ARIMA模型指定参数p、d、q从前文可得d1
选择ARIMA模型的方法 在ACF图中在滞后项为一阶时有一个比较明显的自相关当滞后阶数逐渐增加时偏相关逐渐减小至零。因此可以考虑ARIMA(0,1,1)模型
3拟合模型
使用arima()函数拟合ARIMA模型
arima(ts,orderc(q,d,q))# 拟合ARIMA模型library(forecast)fit-arima(Nile,orderc(0,1,1))fit
# 结果 accuracy(fit)
# 结果 # 指定了d1函数将对时序做一阶差分因此直接将模型应用于原始时序即可。函数返回了移动平均项的系数-0.73以及模型的AIC值。如果还有其他备选模型则可以通过比较AIC值来得到最合理的模型比较的准则是AIC值越小越好。对百分比误差的绝对值做平均的结果是13%。
4模型评估
一般来说一个模型如果合适那么模型的残差应该满足均值为0的正态分布并且对于任意的滞后阶数残差自相关系数都应该为零。也就是说模型的残差应该满足独立正态分布残差间没有关联。
# 模型拟合评估library(ggplot2)## 提取残差df-data.frame(residas.numeric(fit$residuals))## 创建Q-Q图ggplot(df,aes(sampleresid))stat_qq()stat_qq_line()labs(title”Normal Q-Q Plot”)
# 结果 判断序列残差是否满足正态性假设的正态Q-Q图
# 计算值如果服从正态分布则数据中的点会落在图中的线上显然本例的结果还不错## 检验所有滞后阶数的自相关系数是否为零Box.test(fit$residuals,type”Ljung-Box”)
# 结果 # Box.test()函数可以检验残差的自相关系数是否都为零。本例中模型的残差没有通过显著性检验可以认为残差的自相关系数为零。ARIMA模型能够较好地拟合本数据。
5预测
如果模型残差不满足正态性假设或自相关系数为零的假设则需要调整模型、增加参数或改变差分次数。选定模型后就可以用它来做预测。
# 用ARIMA模型做预测forecast()函数forecast(fit,3)
# 结果 autoplot(forecast(fit,3))labs(x”Year”,y”Annual Flow”)
# 结果 用ARIMA(0,1,1)模型对Nile序列做接下来三年的预测。
# 图中的黑线代表点估计浅蓝色和深蓝色区域分别代表80%置信区间和95%的置信区间
6.3 ARIMA模型的自动预测
使用forecast包中的auto.arima()函数实现最优ARIMA模型的自动选取
# ARIMA模型的自动预测使用sunspots数据集library(forecast)fit-auto.arima(sunspots)fit
# 结果 # 函数选定ARIMA模型的参数为p2、d1、q2与其他模型相比这种情况下得到的模型的AIC值最小。
forecast(fit,3)
# 结果 accuracy(fit)
# 结果 # 由于序列中存在值为零的观测值MPE和MAPE这两个准确度度量都失效了这也是这两个统计量的一个缺陷
