后使用快捷导航没有帐号？
用R语言做时间序列分析（附数据集和源码）
查看: 55331|
评论: 0|来自: 数据挖掘DW
摘要: 时间序列（time series）是一系列有序的数据。通常是等时间间隔的采样数据。如果不是等间隔，则一般会标注每个数据点的时间刻度。下面以time series 普遍使用的数据 airline passenger为例。 这是十一年的每月乘客数 ...
时间序列（time series）是一系列有序的数据。通常是等时间间隔的采样数据。如果不是等间隔，则一般会标注每个数据点的时间刻度。下面以time series 普遍使用的数据 airline passenger为例。 这是十一年的每月乘客数量，单位是千人次。如果想尝试其他的数据集，可以访问这里： https://datamarket.com/data/list/?q=provider:tsdl可以很明显的看出，airline passenger的数据是很有规律的。time series data mining 主要包括decompose（分析数据的各个成分，例如趋势，周期性），prediction（预测未来的值），classification（对有序数据序列的feature提取与分类），clustering（相似数列聚类）等。这篇文章主要讨论prediction（forecast，预测）问题。 即已知历史的数据，如何准确预测未来的数据。先从简单的方法说起。给定一个时间序列，要预测下一个的值是多少，最简单的思路是什么呢？（1）mean（平均值）：未来值是历史值的平均。（2）exponential smoothing （指数衰减）：当去平均值得时候，每个历史点的权值可以不一样。最自然的就是越近的点赋予越大的权重。或者，更方便的写法，用变量头上加个尖角表示估计值（3) snaive ： 假设已知数据的周期，那么就用前一个周期对应的时刻作为下一个周期对应时刻的预测值（4） drift：飘移，即用最后一个点的值加上数据的平均趋势介绍完最简单的，下面开始介绍两个time series里面最火的两个强大的算法： Holt-Winters 和 ARIMA。 上面简答的算法都是这两个算法的某种特例。（5）Holt-Winters： 三阶指数平滑Holt-Winters的思想是把数据分解成三个成分：平均水平（level），趋势（trend），周期性（seasonality）。R里面一个简单的函数stl就可以把原始数据进行分解：一阶Holt—Winters假设数据是stationary的（静态分布），即是普通的指数平滑。二阶算法假设数据有一个趋势，这个趋势可以是加性的(additive,线性趋势)，也可以是乘性的(multiplicative,非线性趋势)，只是公式里面一个小小的不同而已。 三阶算法在二阶的假设基础上，多了一个周期性的成分。同样这个周期性成分可以是additive和multiplicative的。 举个例子，如果每个二月的人数都比往年增加1000人，这就是additive；如果每个二月的人数都比往年增加120%，那么就是multiplicative。R里面有Holt-Winters的实现，现在就可以用它来试试效果了。我用前十年的数据去预测最后一年的数据。 性能衡量采用的是RMSE。 当然也可以采用别的metrics：预测结果如下：结果还是很不错的。（6） ARIMA： AutoRegressive Integrated Moving AverageARIMA是两个算法的结合：AR和MA。其公式如下：是白噪声，均值为0， C是常数。 ARIMA的前半部分就是Autoregressive：&，后半部分是moving average：&。 AR实际上就是一个无限脉冲响应滤波器（infinite impulse resopnse）, MA是一个有限脉冲响应（finite impulse resopnse），输入是白噪声。ARIMA里面的I指Integrated（差分）。 ARIMA（p,d,q）就表示p阶AR，d次差分，q阶MA。 为什么要进行差分呢？ ARIMA的前提是数据是stationary的，也就是说统计特性（mean，variance，correlation等）不会随着时间窗口的不同而变化。用数学表示就是联合分布相同：当然很多时候并不符合这个要求，例如这里的airline passenger数据。有很多方式对原始数据进行变换可以使之stationary：（1） 差分，即Integrated。 例如一阶差分是把原数列每一项减去前一项的值。二阶差分是一阶差分基础上再来一次差分。这是最推荐的做法（2）先用某种函数大致拟合原始数据，再用ARIMA处理剩余量。例如，先用一条直线拟合airline passenger的趋势，于是原始数据就变成了每个数据点离这条直线的偏移。再用ARIMA去拟合这些偏移量。（3）对原始数据取log或者开根号。这对variance不是常数的很有效。如何看数据是不是stationary呢？这里就要用到两个很常用的量了： ACF（auto correlation function）和PACF(patial auto correlation function)。对于non-stationary的数据，ACF图不会趋向于0，或者趋向0的速度很慢。 下面是三张ACF图，分别对应原始数据，一阶差分原始数据，去除周期性的一阶差分数据：确保stationary之后，下面就要确定p和q的值了。定这两个值还是要看ACF和PACF：确定好p和q之后，就可以调用R里面的arime函数了。 值得一提的是，R里面有两个很强大的函数： ets 和 auto.arima。 用户什么都不需要做，这两个函数会自动挑选一个最恰当的算法去分析数据。在R中各个算法的效果如下：代码如下：passenger = read.csv('passenger.csv',header=F,sep=' ')p&-unlist(passenger)pt&-ts(p,frequency=12,start=2001)plot(pt)train&-window(pt,start=2001,end=)test&-window(pt,start=2012)library(forecast)pred_meanf&-meanf(train,h=12)rmse(test,pred_meanf$mean) #226.2657pred_naive&-naive(train,h=12)rmse(pred_naive$mean,test)#102.9765pred_snaive&-snaive(train,h=12)rmse(pred_snaive$mean,test)#50.70832pred_rwf&-rwf(train,h=12, drift=T)rmse(pred_rwf$mean,test)#92.66636pred_ses &- ses(train,h=12,initial='simple',alpha=0.2)rmse(pred_ses$mean,test) #89.77035pred_holt&-holt(train,h=12,damped=F,initial="simple",beta=0.65)rmse(pred_holt$mean,test)#76.86677 &without beta=0.65 it would be 84.41239pred_hw&-hw(train,h=12,seasonal='multiplicative')rmse(pred_hw$mean,test)#16.36156fit&-ets(train)accuracy(predict(fit,12),test) #24.390252pred_stlf&-stlf(train)rmse(pred_stlf$mean,test)#22.07215plot(stl(train,s.window="periodic")) &#Seasonal Decomposition of Time Series by Loessfit&-auto.arima(train)accuracy(forecast(fit,h=12),test) #23.538735ma = arima(train, order = c(0, 1, 3), & seasonal=list(order=c(0,1,3), period=12))p&-predict(ma,12)accuracy(p$pred,test) &#18.55567BT = Box.test(ma$residuals, lag=30, type = "Ljung-Box", fitdf=2)
刚表态过的朋友 ()
上一篇：下一篇：
dataguru.cn All Right Reserved.R语言入门（1）时间序列分析
时间序列分析
使用软件为Rstudio，参考，拿手上的数据练习一下时间序列分析。
1、原始数据说明
选择连续9天的数据，共2025条，时间间隔为5分钟。具体情况如下：
2、平稳性检验
所谓平稳，是指因变量围绕着一个常数上下波动。更学术一点，就是是说统计特性（mean，variance，correlation等）不会随着时间窗口的不同而变化。
2.1 时间序列可视化
因为一些高阶函数可能不能准确地产生我想要的图，所以我这里用的是低阶图形命令手动设置label。首先把plot()函数里的默认坐标轴关掉：xaxt='n'，然后通过axis()函数手动添加坐标轴。代码如下：
plot(indata,type = "s",xaxt='n',xlab = "time")
axis(1,c(1,226,451,676,901,1126,1351,1576,1801),c('5:00','5:00','5:00','5:00','5:00','5:00','5:00)
生成图形如下：
由上图可知：
1. 该时间序列存在很强的周期性
2. 该时间序列是平稳的。
还可以用xts函数把原始数据转换成时间序列，也就是把给每个数据对应上一个时间点，代码如下
data_1&-xts(indata_1,seq(ISOdate(0, 0, 0, 5, 0, 0),ISOdate(0, 0, 0, 23, 40, 0),by="5 mins"))
plot(data_1)
进一步绘制每一天的箱型图，代码如下：
library(ggplot2)
ggplot(data,aes(x=factor(STATDATE),y=INDATA)) + geom_boxplot()
由上图可知：
周末的异常点较少，平日的异常点较多，说明平日存在高峰情况。
2.2 自相关和偏自相关图
此外，还可以通过。平稳的序列的自相关图和偏相关图不是拖尾就是截尾。截尾是指在某阶之后，系数都为 0；拖尾是指有一个衰减的趋势，但是不都为0。
固定点上有大于0的值（周期性）
所有点都是0
完全随机（白噪声）
指数衰减到0（拖尾）
正负交替衰减到0（拖尾）
前几个大于0，后面都为0（截尾）
绘制采用R中acf(indata)函数（默认滞后30阶）的自相关图和偏自相关图如下：
由上图可知：
自相关图逐渐衰减到0，存在拖尾性，偏自相关图在5阶后都为0，存在截尾性。说明该时间序列在短期内（30阶）是平稳时间序列。
考虑到滞后阶数为30，而该原始数据的周期性可能要到225阶以后才能体现出来。设置自相关图的滞后阶数为500，代码如下：
acf(indata,lag=500)
pacf(indata,lag=500)
2.3 ADF单位根平稳型检验
单位根检验是指检验序列中是否存在单位根，因为存在单位根就是非平稳时间序列了。[百度百科]
此处使用tseries包中的adf.test()函数对9天数据进行单位根平稳性检验，结果如下：
& adf.test(indata)
Augmented Dickey-Fuller Test
Dickey-Fuller = -7.1486, Lag order = 12, p-value = 0.01
alternative hypothesis: stationary
aternative hypothesis称为备择假设，与原假设相对。这里备择假设是该时间序列平稳。因为p值小于0.05，拒绝原假设，接受备择假设，所以说明该时间序列是平稳的。
有意思的是，我单独拿一天的数据出来做单位根检验，可以发现该时间序列并不平稳。
& adf.test(indata_1)
Augmented Dickey-Fuller Test
Dickey-Fuller = -1.0092, Lag order = 6, p-value = 0.9351
alternative hypothesis: stationary
3、时间序列模型选择
3.1 常用的时间序列模型
下面简单介绍下常用的时间序列模型：
移动平均（MA, moving average）
MA(1): Yt=α0+α1ut-1+ut
MA(q): Yt=α0+α1ut-1+...+αqut-q+ut
自回归模型（AR, Auto regresstion）
AR(1): Yt=β0+β1Yt-1+ut
AR(p): Yt=β0+β1Yt-1+...+βpYt-p+ut
自回归移动平均（ARMA）
ARMA(p, q): Yt=β0+β1Yt-1+...+βpYt-p+α1ut-1+...+αqut-q+ut
差分自回归移动平均（ARIMA）
观测时间序列是否为平稳。若为平稳序列，则用ARMA(p,q)模型；若为非平稳时间序列，则要先进行d阶差分运算后化为平稳时间序列，此处的d即为ARIMA(p,d,q)模型中的d。
3.2 时间序列模型选择
自相关和偏自相关检验：
因为平稳的序列的自相关图和偏相关图不是拖尾就是截尾，所以可根据自相关和偏自相关图的截尾和拖尾来选择时间序列模型。
衰减趋于零(拖尾)
衰减趋于零（拖尾）
ARMA(p, q)
q阶后衰减趋于零（拖尾）
p阶后衰减趋于零（拖尾）
若都拖尾，得到ARMA(p，q)模型，自相关图有几个在两倍标准差之外就能确定p，偏自相关图突出两倍标准差的确定q。可用AIC和SBC准则判断得到的p和q参数值的好坏，这两个指标越小越好。
所幸，R语言的forecast包中提供的auto.arima()函数可以自动找到最优的p，q值，代码如下：
& library(forecast)
& auto.arima(indata,trace = T)
ARIMA(2,0,2) with non-zero mean : 20160.45
ARIMA(0,0,0) with non-zero mean : 24296.42
ARIMA(1,0,0) with non-zero mean : 20279.09
ARIMA(0,0,1) with non-zero mean : 22583.28
ARIMA(0,0,0) with zero mean
: 26869.16
ARIMA(1,0,2) with non-zero mean : 20157.5
ARIMA(1,0,1) with non-zero mean : 20169.87
ARIMA(1,0,3) with non-zero mean : 20159.35
ARIMA(2,0,3) with non-zero mean : 20161.85
ARIMA(1,0,2) with zero mean
ARIMA(0,0,2) with non-zero mean : 21739.44
Best model: ARIMA(1,0,2) with non-zero mean
Series: indata
ARIMA(1,0,2) with non-zero mean
Coefficients:
sigma^2 estimated as 1229:
log likelihood=-10075.54
AIC=20161.08
AICc=20161.11
BIC=20189.15
系数的标准差
sigma^2 estimated
估计值方差
log likelihood
对数似然值
还可以用summary()函数对拟合结果进行描述统计，计算一些误差（Training set error measures:）。
& summary(auto.arima(indata))
Series: indata
ARIMA(1,0,2) with non-zero mean
Coefficients:
sigma^2 estimated as 1229:
log likelihood=-10075.54
AIC=20161.08
AICc=20161.11
BIC=20189.15
Training set error measures: ME
Training set：0. 35.02654 25.4774 -Inf
Inf 0.9615727 -2.
Mean Error
Root Mean Squared Error
Mean Absolute Error
Mean Percentage Error
Mean Absolute Percentage
Mean Absolute Scaled Error
从上述结果可以看出，用来拟合这九天最优的模型是ARIMA(1,0,2) with non-zero mean。with non-zero mean是指常数项不为0，因此我们可以写出该ARIMA模型的公式，如下：
Yt=154.5767+0.9711×Yt-1-0.2664×ut-1-0.0793×ut-2+ut
其中ut、ut-1、ut-2都服从正态分布。
3.3 用模型来预测
在选择好模型后，我们可以使用该模型来做一些预测，例如预测未来一天（225个时段）的数据，代码如下：
fit&-auto.arima(indata)
forecast &- forecast.Arima(fit,h=225,level=c(99.5))
plot.forecast(forecast)
图1 预测未来225个时段（1天）
图2 预测未来7个时段
由图1和图2可知，ARIMA模型仅能适用于短期的预测（7个时段），不能使用于长期的预测（225个时段）。
看过本文的人也看了：
我要留言技术领域：
取消收藏确定要取消收藏吗？
删除图谱提示你保存在该图谱下的知识内容也会被删除，建议你先将内容移到其他图谱中。你确定要删除知识图谱及其内容吗？
删除节点提示无法删除该知识节点，因该节点下仍保存有相关知识内容！
删除节点提示你确定要删除该知识节点吗？用R语言做时间序列分析-土地公文库
用R语言做时间序列分析
用R语言做时间序列分析
代码如下：passenger = #39;#39;,header=F,sep=' ') p&-unlist(passenger)pt&-ts(p,frequency=12,start=2001) plot(pt)train&-window(pt,start=2001,end=)test&-window(pt,start=2012)library(forecast)pred_meanf&-meanf(train,h=12)rmse(test,pred_meanf$mean) #pred_naive&-naive(train,h=12)rmse(pred_naive$mean,test)#pred_snaive&-snaive(train,h=12)rmse(pred_snaive$mean,test)#pred_rwf&-rwf(train,h=12, drift=T)rmse(pred_rwf$mean,test)#pred_ses &- ses(train,h=12,initial='simple',alpha=rmse(pred_ses$mean,test) #第9 / 10页
全文分10页阅读*邮箱不能为空
*密码不能为空
*验证码不能为空
*昵称不能为空
*邮箱不能为空
点击获取验证码
*验证码不能为空
*密码不能为空
*密码不能为空
一个好用的时间序列分析R包astsa
粤ICP备号&&(C)2015版权所有R语言做时间序列分析时,summary给出的结果都是什么意思啊?
R语言做时间序列分析时,summary给出的结果都是什么意思啊?比如下面这个例子：summary(auto.arima(z))Series:z ARIMA(4,0,2) with zero mean Coefficients:ar1 ar2 ar3 ar4 ma1 ma2-0.6 0.5 -0.7s.e.0.8 0.0 0.2sigma^2 estimated as 417.6:log likelihood=-347.56AIC=709.13 AICc=710.73 BIC=725.63Training set error measures:ME RMSE MAE MPE MAPE MASETraining set 2..3 85.53 0.6767153有没有办法检验系数估计的显著性啊?
这个是自动适应参数估计的结果.模型估计为ARIMA(4,0,2),即ARMA(4,2)系数为：ar1 ar2 ar3 ar4 ma1 ma2-0.6 0.5 -0.7s.e.0.8 0.0 0.2s.e.是系数的标准差,系数显著性要自己算,|系数/se| > 1.96 即 95％的置信度sigma^2 estimated 估计值方差log likelihood 对数似然值（这个不用解释了吧）AIC=709.13 AICc=710.73 BIC=725.63再就是下面一堆误差计算ME\x05Mean ErrorRMSE\x05Root Mean Squared ErrorMAE\x05Mean Absolute ErrorMPE\x05Mean Percentage ErrorMAPE\x05Mean Absolute PercentageMASE\x05Mean Absolute Scaled Erro
我有更好的回答：
剩余:2000字
与《R语言做时间序列分析时,summary给出的结果都是什么意思啊?》相关的作业问题
egression,有本‘例解回归分析’真心不错,可惜和R语言没什么关系；用R专门讲回归分析的真心不多（至少我没见到过多少）,但随便找本R的教材,里面都至少会有一章讲到regression. 再问： 真嘚咧？
现在完成时 in是一段时间
首先,Matlab为收费软件价格不低,安装文件比较大.但是随着发展,以及R语言等免费软件的兴起,Matlab价格出现多个版本商业、教育、学生、个人版.学生版功能与商业版一致,但价格非常优惠了不到100美元.R语言为开源免费软件,安装文件比较小,而且有着很多牛人为其提供工具包,发展迅速.在图像展现方面似乎以及高于Matl
She is getting better day by day.她的健康状况一天天好起来用进行时
day and night, day by day, day after day等表示时间渐进的词作状语时,谓语动词可以用一般现在时、现在进行时或一般过去时.
EVIEWS比较好,SPSS现在都不是STATA的对手了
数据在哪里需要文字分析吗我的邮箱发给你了把数据和借鉴的论文一起发来哈
spss ---分析----预测 ----创建模型 自变量不用管,将你的数量变量 移入因变量中 方法使用专家建模器保存 里面有个 预测值 选中 点确定 就可以出来创建的时间序列了
做时间序列分析,最强大最方便的是EViews,包括单位根检验、VAR模型、协整检验等等.需要的话,数据发给我,我可以帮您.
什么叫小值波动?利用eviews分析时间序列,1首先应绘制时序图,分析波动和趋势2做相关系数图,3做单位根检验,根据相关系数判断数据平稳性如果数据序列平稳,则根据相关系数图判断模型形式和阶数不平稳的话则需要差分等方法将数据进行平稳化处理平稳化后的序列应该通过单位根检验,再根据相关系数图判断模型形式和阶数
在地球表面的物体受到的重力等于万有引力为：mg=GMmR2，得：GM=R2g根据卫星受到的万有引力等于向心力得：GMm(2R)2=mω2o(2R)，得：ω=GM8R3=g8R=12g2R又GMm(2R)2=ma，得：a=GM4R2=14g故答案为：12g2R，14g．
大部分的时间序列数据都是非平稳的；自相关应该是一个前提假设吧
自相关系数拖尾,偏自相关截尾 做AR模型自相关系数截尾,偏自相关拖尾 做MA模型你这个要做ARIMA模型了,至于滞后期数,要一个一个试验 找AIC SC 最小的期数
你怎么又做欧氏距离,有最小距离法?啥意思啊,你是要聚类吗? 再问： 嗯 聚类 首先得把这些数据的欧氏距离做出来 再答： 你是不知道怎么读文件？还是不知道怎么用R来做欧式距离？用R做距离的话用dist函数就行了。再问： 读取用的是read.table 公式是dist(rbind(A,B),method="euclidea
时间序列建模器 图表那个选项卡 左下勾选 拟合值 就可以了.我的为什么不出现预测值啊啊啊啊~
第二个：设定17,使每次结果一样,去掉响应变量species,用于做模型,最后画出图 再答： 通过？randomforest使用帮助文件再问： 它那里写的是“The `unsupervised' case”，请问是什么意思？ 再答： 非监督案例。算法分为：监督和非监督类型，聚类属于非监督再问： Do MDS on 1
1平时看文字注意训练速度2做试卷时不要太犹豫,你多犹豫几分钟的答案和不犹豫的答案相差不远3写字尽量写快点,可以请教写得又快又好的同学4整张试卷要有一个时间的大体分配,最后作文部分一定留一个小时到40分钟,这是之前的（根据你写作文的速度吧.因人而异,我就试过30分钟不到做完作文的）
在电场中运动那部分时间设为 t电,可认为粒子是从速度为0开始被电场加速到V而磁场对粒子只起偏转作用,不影响粒子的速度大小.所以　qU＝mV^2 / 2(qE)*t电＝mVE＝U / d得　t电＝2*d / V电场每加速一次粒子,就可使粒子获得动能 q*U所以共加速的次数是　N＝（m V^2 / 2）/ (qU)＝m V
楼主的困惑源自两个方面：1.混淆了研究对象；2.没能深刻体会到加速度是一个比值,单纯考虑力或质量的大小都明确不了加速度的大小.加速度的公式是a=F/m,很明显,这是一个比值.一台汽车和一台自行车做自由落体时的加速度相等,楼主受困惑的是,明明这两个东西所受到的引力不一样啊?为什么还会相等呢?没错,他们受到的引力是不一样,}