如何用r将时间序列分析ppt解成趋势,季节性和误差成分

用神经网络时间序列做预测,预测结果整体还好基本误差都非常小,但是偶尔的几个预测值的误差大的也离谱_百度知道
用神经网络时间序列做预测,预测结果整体还好基本误差都非常小,但是偶尔的几个预测值的误差大的也离谱
误差如:【0.03 0.06 0. 0.03 0.19 0.】为什么会这样啊
我有更好的答案
0.0相差1.5W倍,这肯定是不行的,太不稳定。看看是不是忘记对数据进行归一化?没归一化的话,会导致数量级大的输入的权值占主导地位,弱化其他输入向量维的作用。如果不是归一化的原因,看看是不是网络结构有问题,例如改变隐层节点数、改变输入向量结构,或者干脆换种神经网络。
我是用 matlab2012b自带的网络工具箱做的,只输入了数据然后经过多次调解训练数据测试数据,还有隐藏网络数节点数之后自动生成的,有什么缺陷吗
你归一化了吗?尽量少用nntool,直接用newff函数这样的代码来建立网络较好,便于控制参数。
采纳率:77%
为您推荐:
其他类似问题
时间序列的相关知识
换一换
回答问题,赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。用数据预测未来:时间序列分析 -ZAKER新闻
人人都是产品经理
对于本文内容,小编只知道作者介绍了一种用数据预测未来的方法——时间序列分析。……嗯,内容灰常灰常灰常烧脑,各位看官 enjoy~应用背景:通过分析序列进行合理预测,做到提前掌握未来的发展趋势,为业务决策提供依据,这也是决策科学化的前提。时间序列分析:时间序列就是按时间顺序排列的一组数据序列。时间序列分析就是发现这组数据的变动规律并用于预测的统计技术。分析工具:SPSS(数据分析的重量级应用,与 SAS 二选一)实践案例:通过历史数据预测未来数据,所涉及的都是最简单的实践,抛砖引玉,重在方法,不论多复杂的数据,方法是一样的。如已知前几年每月的销售量,预测未来的销售量。一、时间序列分析简介时间序列分析有三个基本特点:假设事物发展趋势会延伸到未来预测所依据的数据具有不规则性不考虑事物发展之间的因果关系并不是所有的时间序列都一定包含四种因素,如以年为单位的诗句就可能不包含季节变动因素。四种因素通常有两种组合方式:四种因素相互独立,即时间序列是四种因素直接叠加而成的,可用加法模型表示:
Y=T+S+C+I四种因素相互影响。即时间序列是四种因素相互综合的结果,可用乘法模型表示:Y=T*S*C*I其中,原始时间序列值和长期趋势可用绝对数表示;季节变动、循环变动、不规则变动可用相对数(变动百分比)表示。二、季节分解法当我们对一个时间序列进行预测时,应该考虑将上述四种因素从时间序列中分解出来。为什么要分解这四种因素?分解之后,能够克服其他因素的影响,仅仅考量一种因素对时间序列的影响。分解之后,也可以分析他们之间的相互作用,以及他们对时间序列的综合影响。当去掉这些因素后,就可以更好的进行时间序列之间的比较,从而更加客观的反映事物变化发展规律。分解之后,序列可以用来建立回归模型,从而提高预测精度。所有的时间序列都要分解这四种因素吗?通常情况下,我们考虑进行季节因素的分解,也就是将季节变动因素从原时间序列中去除,并生成由剩余三种因素构成的序列来满足后续分析需求。为什么只进行季节因素的分解?时间序列中的长期趋势反映了事物发展规律,是重点研究的对象;循环变动由于周期长,可以看做是长期趋势的反映;不规则变动由于不容易测量,通常也不单独分析。季节变动有时会让预测模型误判其为不规则变动,从而降低模型的预测精度综上所述:当一个时间序列具有季节变动特征时,在预测值钱会先将季节因素进行分解。步骤:定义日期标示变量:即先将序列的时间定义好,才能分析其时间特征。了解序列发展趋势:即序列图,确定乘性还是加性进行季节因素分解建模分析结果解读预测1、定义日期标示变量时间序列的特点就是数据根据时间点的顺序进行排列,因此分析之前,SPSS 需要知道序列的时间定义,然后才能进行分析时间特征。根据源数据的格式进行选择,并输入第一个个案的具体数值。此时会在源文件中生成三个新的变量。2、了解序列发展趋势完成日期标示变量的定义之后,需要先对时间序列的变化趋势有所了解,便于选择合适的模型。即通过序列图,确定模型是乘性还是加性。变量为 " 销售数据 ",时间轴标签为 "DATE – ",也就是我们自定义的时间。数据销量序列图如何根据序列图来判断模型的乘性或加性?如果随着时间的推移,序列的季节波动变得越来越大,则建议使用乘法模型。如果序列的季节波动能够基本维持恒定,则建议使用加法模型。本例很明显:随着时间变化,销售数据的季节波动越来越大,那么使用乘法模型会更精确。3、进行季节因素分解变量为 " 销售数据 ",且根据序列图我们知道时间序列模型为乘性。提示您会新生成四个变量ERR(误差序列):从时间序列中移除季节因素、长期趋势、和循环变动之后留下的序列,也就是原始序列中的不规则变动构成的序列。SAS(季节因素校正后序列):是移除原始序列中的季节因素后的校正序列。SAF(季节因子):是从序列中分解出的季节因素。其中的变量值根据季节周期的变动进行重复,如本例中季节周期为 12 个月,所以这些季节因子没 12 个月重复一次。STC(长期趋势和循环变动趋势):这是原始序列中长期趋势和循环变动构成的序列。如图,周期为 12 个月,季节因子 12 个月循环一次。完成季节因素分解后的序列和原始序列之间有什么差异?通过回执序列图的方法把原始序列和除去季节因子的三个序列(误差序列、季节因素校正后序列、长期无视和循环变动序列)进行比较。要做四个序列图,会有四个变量:原始序列:使用变量 " 销售数据 ";误差序列:使用变量 "ERR";季节因素校场后序列:使用变量 "SAS"长期趋势和循环变动序列:使用变量 "STC"蓝色线:原始序列紫色线:长期趋势和循环变动序列浅棕色:季节因素校正后序列绿色线:误差序列(不规则变动)因为误差序列数值非常小,所以长期趋势和循环变动序列(长期趋势 + 循环变动)与季节因素校正后序列(长期趋势 + 循环变动 + 不规则变动,即误差)能够基本重合。在单独做 " 季节因子 SAF" 的序列图:因为是做 " 季节因子 " 的序列图,所以只有一个变量 " 季节因子 SAF"我们看出:季节因素的周期是 12 个月,先下降,然后上升到第一个顶点,再有略微的下降后,出现明显的上升趋势,到第七个月时达到峰值,然后一路下跌,直到最后一个月份有所回升,之后进入第二个循环周期。通过对原始序列的季节分解,我们更好的掌握了原始序列所包含的时间特征,从而选用适当的模型进行预测。三、专家建模法时间序列的预测步骤有四步:绘制时间序列图观察趋势分析序列平稳性并进行平稳化时间序列建模分析模型评估与预测平稳性主要是指时间序列的所有统计性质都不会随着时间的推移而发生变化。对于一个平稳的时间序列,具备以下特征:均数和方差不随时间变化自相关系数只与时间间隔有关,与所处的时间无关自相关系数是研究序列中不同时期的相关系数,也就是对时间序列计算其当前和不同滞后期的一系列相关系数。平稳化的方法——差分。差分就是指序列中相邻的两期数据之差。一次差分 =Yt-Yt-1二次差分 = ( Yt-Yt-1 ) - ( Yt-1-Yt-2 ) 具体的平稳化操作过程会有专家建模法自动处理,我们只需要哼根据模型结果独处序列经过了几阶差分即可。时间序列分析操作:要分析所有变量,所以选择 " 销售数据 "。【专家建模器】–【条件】,勾选 " 专家建模器考虑季节性模型 "。勾选 " 预测值 ",目的是生成预测值,并保存模型。时间序列分析结果解读该表显示了经过分析得到的最优时间序列模型及其参数,最优时间 U 型猎魔性为 ARIMA(0,1,1)(0,1,1)求和自回归移动平均模型 ARIMA(p,d,q) ( P,D,Q ) p:出去季节性变化之后的序列所滞后的 p 期,通常为 0 或 1,大于 1 的情况很少;d:除去季节性变化之后的序列进行了 d 阶差分,通常取值为 0,1 或 2;q:除去季节性变化之后的序列进行了 q 次移动平均,通常取值 0 或 1,很少会超过 2;P,D,Q 分别表示包含季节性变化的序列所做的事情。因此本例可解读为:对除去季节性变化的序列和包含季节性变化的序列分别进行了一阶差分和一次移动平均,综合两个模型而建立出来的时间序列模型。该表主要通过 R 方或平稳 R 方来评估模型拟合度,以及在多个模型时,通过比较统计量找到最优模型。由于原始变量具有季节性变动因素,所以平稳的 R 方更具有参考意义,等于 32.1%,拟合效果一般。该表提供了更多的统计量可以用来评估时间序列模型的拟合效果。虽然平稳 R 方仅仅是 32.1%,但是 " 杨 - 博克斯 Q(18)" 统计量的显著性 P=0.706,大于 0.05(此处 P&0.05 是期望得到的结果),所以接受原假设,认为这个序列的残差符合随机分布,同时没有离群值出现,也都反映出数据的拟合效果还可以接受。时间序列应用预测:未来一年是到 2016 年 12 月,手动输入即可。这是未来一年的销售趋势。如果想从全局来观察预测趋势,可以在把这一年的趋势和以前的数据连接起来此时的变量应该是 " 原始的销售数量 " 和 "2016 年的预测销售数量 "。结果如下:也可以在表中查看具体的数值:作者:膝盖哥,是一枚 " 跪着提需求 " 的产品经理。常说 " 不用不用,真的不用了,我跪着就好!"本文由 @膝盖哥 原创发布于人人都是产品经理。未经许可,禁止转载。
相关标签:
原网页已经由 ZAKER 转码排版
猎云网14小时前
猎云网15小时前
36氪6小时前
36氪7小时前
36氪14小时前
36氪9小时前
猎云网19小时前
亿欧网4小时前
知乎每日精选19小时前
砍柴网6小时前
Tech2IPO6小时前
砍柴网6小时前
砍柴网6小时前
砍柴网6小时前
36氪11小时前 上传我的文档
 上传文档
 下载
 收藏
粉丝量:115
该文档贡献者很忙,什么也没留下。
 下载此文档
第05章 多元时间序列分析方法
下载积分:30
内容提示:第05章 多元时间序列分析方法
文档格式:PDF|
浏览次数:306|
上传日期: 02:07:20|
文档星级:
全文阅读已结束,如果下载本文需要使用
 30 积分
下载此文档
该用户还上传了这些文档
第05章 多元时间序列分析方法
关注微信公众号23分享收藏文章被以下专栏收录Time series Analysis with R
原文例子是来自:
&&时间序列分析及应用&&
&&R语言时间序列中文教程&&
&&A Little Book of R for Time Series&&&
&&Analysis of Integrated and Cointegrated Time Series with R&&
1. examples
1.1 洛杉矶年降水
install.packages(&TSA&)
library(TSA)
win.graph(width = 4.875, height = 2.5, pointsize = 8) #产生一个win界面
data(larain)
plot(larain, ylab = 'Inches', xlab = 'Year', type = 'o')
win.graph(width = 3, height = 3, pointsize = 8)
plot(y = larain, x = zlag(larain), ylab = 'Inches', xlab = 'Previous Year Inches')
1.2 &化工过程
win.graph(width = 4.875, height = 2.5, pointsize = 8)
data(color)
plot(color, ylab = 'Color Property', xlab = 'Batch', type = 'o')
win.graph(width = 3, height = 3, pointsize = 8)
plot(y = color, x = zlag(color), ylab = 'Color Property', xlab = 'Previous Batch Color Property') #zlag函数将对象后移n位
1.3 加拿大野兔年度丰度
win.graph(width = 4.875, height = 2.5, pointsize = 8)
data(hare)
plot(hare, ylab = 'Abundance', xlab = 'Year', type = 'o')
1.4 建模策略
主要步骤:
(1)模型识别
(2)模型拟合
(3)模型诊断
2. 基本概念
2.1 时间序列与随机过程
2.2 均值、方差和协方差
随机余弦波
3.1 确定性趋势与随机趋势
随机趋势:对完全一样的过程进行多次模拟,可能展现出完全不同的”趋势“,即为随机趋势。
确定趋势:在一事实上的时间周期内,可以预测某一过程的发生情况的趋势,即为确定趋势。
3.2 常数均值的估计
3.3 回归方法
时间的线性趋势
时间的二次趋势
周期性或季节性趋势
data(tempdub)
head(tempdub)
month. = season(tempdub)
head(month.)
model2 = lm(tempdub ~ month. - 1) # -1删除intercept term
summary(model2)
model3 = lm(tempdub ~ month.)
summary(model3)
har. = harmonic(tempdub, 1)
model4 = lm(tempdub ~ har.)
summary(model4)
win.graph(width = 4.875, height = 2.5, pointsize = 8)
plot(ts(fitted(model4), freq = 12, start = c(1964, 1)), ylab = 'Temperature', type = 'l', ylim = range(c(fitted(model4), tempdub)))
points(tempdub)
3.4 回归估计的可靠性和有效性
最小二乘估计量
最佳线性无偏估计量(BLUE)
广义最小二乘估计量(GLS)
3.5 回归结果的解释
残差的标准差
决定系数R^2 : 观测序列与估计的趋势之间的样本相关系数的平方,也表示序列的变化被估计的趋势所解释的部分。
3.6 残差分析
无法观测的随机项{Xt}可以通过残差来估计或预测
plot(y = rstudent(model3), x = as.vector(time(tempdub)),
xlab = 'Time', ylab = 'Standardized Residuals', type = 'o')
plot(y = rstudent(model3), x = as.vector(fitted(model3)), xlab = 'Fitted Trend Values',
ylab = 'Standarized Residuals', type = 'n')
points(y = rstudent(model3), x = as.vector(fitted(model3)), pch = as.vector(season(tempdub)))
hist(rstudent(model3), xlab = 'Standarized Residual')
win.graph(width = 2.5, height = 2.5, pointsize = 8)
qqnorm(rstudent(model3))
QQ图:显示了数据的贫僧数和根据正态分布计算的理论分位数;对于正态分布数据,QQ图看起来近似于一条直线。
Shapiro-Wilk检验,本质是计算残差与其相应的正态分位数之间的相关系数,相关系数越小,就越有理由否定正态性。
样本自相关函数&
#季节均值模型残差的样本自相关系数
win.graph(width = 4.875, height = 3, pointsize = 8)
acf(rstudent(model3))
4. 平稳时间序列模型
4.1 一般线性过程
一般线性过程{Yt}可表示成现在和过去白噪声变量的加权线性组合:
4.2 滑动平均过程
当有限个系数不为零的时候,得到所谓的滑动平均过程
q阶滑动平均过程(MA(q)): 过程大于1阶滞后时,不存在自相关
一阶滑动平均过程:不存在自相关
二阶滑动平均过程:
一般MA(q)过程:
4.3 自回归过程
自回归过程是用自身做回归变量
p阶自回归过程{Yt}满足方程:
一阶自回归过程:
AR(1)模型的二阶滞后也具有强正自相关
AR(1)过程的平稳性:
& 通常称为AR(1)过程的平稳条件
二阶自回归过程:
AR(2)过程的平稳性:
AR(2)过程的自相关函数:
Yule-Walker方程
一般自回归过程:
4.4 自回归滑动平均混合模型
{Yt}称为自回归滑动平均混合过程,阶数分别为p和q
4.5 可逆性
特征多项式
win.graph(width = 4, height = 3, pointsize = 8)
data(ma2.s)
plot(ma2.s, ylab = expression(Y[t-1]), type = 'o')
&span style=&color:#ff6666;&&# Y[t-1]表示下标为t-1的Y&/span&
5. 非平稳时间序列模型
具有时变均值的任何时间序列都是非平稳的,形如:
5.1 通过差分平稳化
plot(diff(log(oil.price)), main = '石油价格序列取对数后的差分时序图', ylab = 'Change in Log(Price)', type = 'l')
5.2 ARIMA模型
若一个时间序列{Yt}的d次差分是一个平衡的ARMA过程,则称{Yt}为自回归滑动平均求和模型
若Wt服从ARMA(p,q)模型,我们称{Yt}是ARIMA(p,d,q)过程
IMA(1,1)模型:
IMA(2,2)模型:&& & &即&
data(ima22.s)
plot(ima22.s, main = 'IMA(2,2)',ylab = 'IMA(2,2) Simulation', type = 'o')
plot(diff(ima22.s), main = 'IMA(2,2)序列的一次差分', ylab = 'First Difference', type = 'o')
plot(diff(ima22.s, main = 'IMA(2,2)序列的二次差分', difference = 2), ylab = 'Differenced Twice', type ='o')
ARI(1,1)模型:
5.4 其他变换
百分比变动和对数
幂变换:只有取正值的数据才能使用幂变换
Box-Cox变换
6.模型识别
6.1 样本自相关函数的性质
对于MA(q)模型,当滞后超过q时,自相关函数为0。因此,样本自相关函数是MA过程阶数的一个良好指示器,但AR(p)模型的自相关函数在琔滞后之后不会变为0,它们是逐渐衰减而不是突然截断。
6.2 &偏自相关函数和扩展的自相关函数
偏自相关函数被定义为预测误差之间的相关系数
样本偏自相关函数(PACF)
样本ACF和PACF为识别纯AR(p)或MA(q)模型提供了有效的工具,但对混合 ARMA模型来说,其理论ACF和PACF有着无限多的非零值,使得根据样本ACF和PACF来识别模型非常困难。
ACF指自动关系性,ACF即Auto-CorrelationFunction的简称. 比方说,股票价格今天的价格跟昨天的价格有关系,明天的价格会跟今天的或者昨天的价格有关系。它们之间的关系性便用ACF来衡量。
PACF被称作不完全自动关系性。自动关系性ACF中存在着线性关系性和非线性关系性。不完全自动关系性就是把线性关系性从自动关系性中消除。如果在线性关系性被去除以后,两个时间点之间的关系性也就是不完全关系性。当PACF近似于0,这表明两个时间点之间的关系性是完全由线性关系性所造成的。如果不完全关系性在两个时间点之间不近似于0,这表明线性模型是不能够表达这两个时间点之间的关系。
6.3 非平稳性
Dickey-Fuller单位根检验
ARMA模型ACF和PACF的一般特征
ARMA(p,q) , p &0, q&0
滞后q阶后截尾
滞后p阶后截尾
7. 参数估计
7.1 矩估计
通过令样本矩等于相应的理论矩,并求解所得方程以求得任意未知参数的估计。
滑动平均模型
噪声方差估计
7.2 最小二乘估计
自回归模型
滑动平均模型
7.3 极大似然与无条件最小二乘
极大似然估计
对数似然函数
无条件最小二乘
7.4 估计的性质
8. 模型诊断
8.1 残差分析
残差的正态性
残差的自相关
Ljung-Box检验
8.2 过度拟合和参数冗余&
10. 季节模型&
10.1 季节ARIMA模型
16.R语言时间序列中文教程实例
(1)plot、plot.ts、ts.plot
# sequence of integers from -5 to 5
y = 5*cos(x)
par(mfrow=c(3,2))
# multifigure setup: 3 rows, 2 cols
plot(x, main=&plot(x)&)
# 如果数据是时间序列对象,使用plot()命令就足够了
plot(x, y, main=&plot(x,y)&)
#如果数据是平常序列,使用plot.ts()也可以做时间绘图
plot.ts(x, main=&plot.ts(x)&)
plot.ts(x, y, main=&plot.ts(x,y)&)
ts.plot(x, main=&ts.plot(x)&)
ts.plot(ts(x), ts(y), col=1:2, main=&ts.plot(x,y)&)
# note- x and y are ts objects
(2)正态性检验、lag.plot、stl结构拆分
dljj = diff(log(jj))
#log及差分处理&div style=&background:#EAEAEA;&&&p align=&left&&#正态性检验&/p&&p align=&left&&shapiro.test(dljj)          #&span style=&color:#005555;&& test for normality&/span& 测试结果的正态分布的性质&/p&&p align=&left&&&div style=&background:#EAEAEA;&&&p align=&left&&par(mfrow=c(2,1))        #&span style=&color:#005555;&& set up the graphics &/span&&span style=&color:#005555;&&设置为两图的输出&/span&&/p&&p align=&left&&hist(dljj, prob=TRUE, 12)   #&span style=&color:#005555;&& histogram&/span&柱形分布图&span style=&color:#005555;&&    &/span&&/p&&p align=&left&&lines(density(dljj))     #&span style=&color:#005555;&& smooth it - ?density for details&/span&柱形分布图的曲线 &span style=&color:#005555;&& &/span&&/p&&p align=&left&&qqnorm(dljj)             #&span style=&color:#005555;&& normal Q-Q plot  QQ&/span&&span style=&color:#005555;&&图&/span&&/p&&p align=&left&&qqline(dljj)             #&span style=&color:#005555;&& add a line    &/span&&span style=&color:#005555;&&在&/span&&span style=&color:#005555;&& QQ&/span&&span style=&color:#005555;&&图&/span&&span style=&color:#005555;&&上&/span&&span style=&color:#005555;&&加&/span&&span style=&color:#005555;&&直线&/span&&/p&&/div&
&/p&&/div&延迟图表,也就是lag.plot
lag.plot(dljj, 9, do.lines=FALSE)  #do.lines是否画线
&/pre&&pre code_snippet_id=&485640& snippet_file_name=&blog__6929768& name=&code& class=&plain&&&span style=&color:#005555;&&&/span&&div style=&background:#EAEAEA;&&&p align=&left&&#ACF与PACF&/p&&p align=&left&&par(mfrow=c(2,1)) #&em&&span style=&color:#005555;&& The power of accurate observation is commonly called cynicism&/span&&/em& &/p&&p align=&left&&                  #&em&&span style=&color:#005555;&&     by those whohave not got it.&/span&&/em&&span style=&color:#005555;&& - George Bernard Shaw&/span&&/p&&p align=&left&&acf(dljj, 20)     #&span style=&color:#005555;&& ACF to lag 20 - no graph shown... keep reading&/span&&/p&&p align=&left&&pacf(dljj, 20)    #&span style=&color:#005555;&& PACF to lag 20 - no graph shown... keep reading&/span&&/p&&p align=&left&&#&span style=&color:#005555;&& !!NOTE!! acf2 onthe line below is NOT available in R... details follow the graph below&/span&&/p&&p align=&left&&acf2(dljj)        #&span style=&color:#005555;&& this is what you'll see below&/span&&/p&&/div&
&p align=&left&&#stl结构拆分&/p&&p align=&left&&Log(jj)=趋势+季节+误差&/p&&p align=&left&& log(jj) =trend + season + error&/p&结构拆析的R命令中stl(), 下面语句中stl命令中输入的是lag变型后的jj数据。
&/pre&&p&&/p&&p&&img src=&https://img-blog.csdn.net/55638?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvbGlsYW5mZW5nMTk5MQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast& alt=&& /&&/p&&p&&/p&&p&&/p&&p&&/p&&pre code_snippet_id=&485640& snippet_file_name=&blog__1092693& name=&code& class=&plain&& #生成模型
Q = factor(rep(1:4,21))
# make (Q)uarter factors [that's repeat 1,2,3,4, 21 times]
trend = time(jj)-1970
# not necessary to &center& time, but the results look nicer
reg = lm(log(jj)~0+trend+Q, na.action=NULL)
# run the regression without an intercept
#-- the na.action statement is to retain time series attributes
summary(reg)
lm(formula = log(jj) ~ 0 + trend + Q, na.action = NULL)
Residuals:
-0.262 -0.060
Coefficients:
Estimate Std. Error t value Pr(&|t|)
trend 0.167172
&2e-16 ***
&2e-16 ***
&2e-16 ***
&2e-16 ***
&2e-16 ***
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1254 on 79 degrees of freedom
Multiple R-squared: 0.9935,
Adjusted R-squared: 0.9931
F-statistic:
2407 on 5 and 79 DF,
p-value: & 2.2e-16
&/pre&&pre code_snippet_id=&485640& snippet_file_name=&blog__206088& name=&code& class=&plain&&##&span style=&font-family: Arial, Helvetica, sans-&&上面语句第一行用来生成所需要的Q函数。第二行用来生成线性回归模型中的时间输入,然后存储于叫做trend的数列中。第三行语句是用lm()命令来建立线性模型。在R语言中,~用来分隔进行预测的变量和被预测的变量。&/span&&/pre&&pre code_snippet_id=&485640& snippet_file_name=&blog__1021726& name=&code& class=&plain&&#接下来要进行的工作是,将预测的数据和实观察的数据进行比较
&div style=&background:#EAEAEA;&&&p align=&left&& plot(log(jj),type=&o&)   #&span style=&color:#005555;&& the data in black with little dots&/span& &/p&&p align=&left&&    lines(fitted(reg),col=2) #&span style=&color:#005555;&& the fitted values in bloody red - or use &/span&lines(reg$fitted, col=2)&/p&&/div&
#&span style=&font-size: 12 font-family: 宋体;&&用误差的绘图来确定误差的变化是否比较小范围之内&/span&&div style=&background:#EAEAEA;&&&p align=&left&&  par(mfrow=c(2,1))&/p&&p align=&left&&    plot(resid(reg))      #&span style=&color:#005555;&& residuals - &/span&reg$resid&span style=&color:#005555;&& is same as &/span&resid(reg) &/p&&p align=&left&&     acf(resid(reg),20)    #&span style=&color:#005555;&& acf of the resids&/span& &/p&&/div&
ts.intersect
combine 数据ts.intersect()命令把mort 数据、part数据以及part经过延迟四周的数据捆绑在了一起&div style=&background:#EAEAEA;&&&p align=&left&& ded =ts.intersect(mort,part,part4=lag(part,-4),  dframe=TRUE) #&span style=&color:#005555;&& tie them together in a data frame&/span&&/p&&/div&&div style=&background:#EAEAEA;&&&p align=&left&&  fit =lm(mort~part+part4, data=ded, na.action=NULL)         #&span style=&color:#005555;&& now the regression will work&/span&&/p&&/div&&div style=&background:#EAEAEA;&&&p align=&left&&&span style=&color:#005555;&& &/span&&span style=&font-family: Arial, Helvetica, sans- background-color: rgb(240, 240, 240);&&summary(fit)  #&span style=&font-family: white-space: color: rgb(0, 85, 85); background-color: rgb(234, 234, 234);&&接下来使用&/span&&span style=&font-family: white-space: color: rgb(0, 85, 85); background-color: rgb(234, 234, 234);&&summary&/span&&span style=&font-family: white-space: color: rgb(0, 85, 85); background-color: rgb(234, 234, 234);&&命令输出所有估计完成的参数。&/span&&/span&&/p&&p align=&left&&&span style=&font-family: Arial, Helvetica, sans- background-color: rgb(240, 240, 240);&&&/span&&span style=&font-family: Arial, Helvetica, sans- background-color: rgb(240, 240, 240);&&  &/span&&/p&&/div&
(3)ARIMA(自回归移动平均集成模型)
ARIMA模型。AR指自动回归,MA指移动平均,I指集成。
&pre name=&code& class=&plain&&&strong style=&background-color: rgb(255, 204, 204);&&arima.sim&/strong&
x1 &- arima.sim(list(order = c(1, 0, 0), ar = 0.9), n = 100)
x2 &- arima.sim(list(order = c(1, 0, 0), ar = -0.9), n = 100)
par(mfrow = c(1, 2))
&strong style=&background-color: rgb(255, 204, 204);&&acf与pacf&/strong&par(mfcol=c(2,2))
acf(x1, 20)
acf(x2, 20)
pacf(x1, 20)
pacf(x2, 20)&p align=&left&&&strong&分析:&/strong&&/p&&p align=&left&&&strong&对数列x1来讲,ACF从lag1到lag12都高于蓝色虚线,也就是说两个时间点的距离在1到12之间它们的自动关系都是正面的。所有线性关系在x1数列中被清除了结果是Partial ACF的x1图表。在PACFx1图表上可以看到只有lag0值为1,其他的lag上的关系值都低于蓝色虚线,近似于0,也就是说在x1数列中存在的自动关系基本上是线性关系。使用线性模型可以符合x1数列的要求。&/strong&&/p&&p align=&left&&&strong& &/strong&&/p&&p align=&left&&&strong&对数列x2来讲,两个时间点之间的ACF关系有负面的也有正面的,但经过去除线性关系以后所有不完全自动关系都接近于0,也就是说x2数列中的自动关系也基本上是线性的。使用线性模型符合x2数列要求。&/strong&&/p&
&strong style=&background-color: rgb(255, 204, 204);&&移动平均MA模型&/strong&&span style=&background-color: rgb(255, 204, 204);&&&/span&&pre name=&code& class=&plain& style=&font-weight:&&x = arima.sim(list(order=c(0,0,1), ma=.8), n=100)
par(mfcol=c(3,1))
plot(x, main=(expression(MA(1)~~~theta==.8)))
pacf(x,20)分析:&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&MA1&/span&&span style=&font-size: 12 font-family: 宋体;&&代表了数列输出的结果。最上面显示的是数列的数目,横向坐标从&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&1&/span&&span style=&font-size: 12 font-family: 宋体;&&到&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&100&/span&&span style=&font-size: 12 font-family: 宋体;&&代表有&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&100&/span&&span style=&font-size: 12 font-family: 宋体;&&个数据。中间的图显示的是&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&x&/span&&span style=&font-size: 12 font-family: 宋体;&&数列的自动关系性,在&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&lag0&/span&&span style=&font-size: 12 font-family: 宋体;&&上有很强的自动关系性,因为&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&Lag0&/span&&span style=&font-size: 12 font-family: 宋体;&&表示没时间延迟,当前数据与它有完全的关系性。所以&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&lag0&/span&&span style=&font-size: 12 font-family: 宋体;&&上的关系性什么都不代表。而&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&Lag1&/span&&span style=&font-size: 12 font-family: 宋体;&&上存在着正面的关系性,即当前数据与前一个数据的关系性是正面的。其他延迟上的自动关系性都近似于&/span&&span lang=&EN-US& style=&font-size: 12 font-family: 'Courier New';&&0&/span&&span style=&font-size: 12 font-family: 宋体;&&,可以不被考虑。&/span&&strong&
&span style=&background-color: rgb(255, 204, 204);&&&strong&# an AR2 &/strong&&/span&&span style=&background-color: rgb(255, 204, 204);&&&strong&x = arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=100) 
par(mfcol=c(3,1))
plot(x, main=('expression(AR(2)~~~phi[1]==1~~~phi[2]==-.9)'))
acf(x, 20)
pacf(x, 20)
&/strong&&/span&&span style=&background-color: rgb(255, 204, 204);&&&strong&
&/strong&&/span&
&strong style=&background-color: rgb(255, 204, 204);&&# an ARIMA(1,1,1) &/strong&&strong style=&background-color: rgb(255, 204, 204);&&&/strong&&pre name=&code& class=&plain&&x = arima.sim(list(order=c(1,1,1), ar=.9, ma=-.5), n=200)
par(mfcol=c(3,1))
plot(x, main=('expression(ARIMA(1,1,1)~~~phi==.9~~~theta==-.5)'))
acf(x, 30)
# the process is not stationary, so there is no population [P]ACF ...
pacf(x, 30) # but look at the sample values to see how they differ from the examples above
&strong&自动回归移动平均集成模型的估计分析,即ARIMA模型&/strong&&strong&x = arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=100) # simulate some data
(x.fit = arima(x, order = c(1, 0, 1)))   # fit the model and print the results
&/strong&&strong&Call:
arima(x = x, order = c(1, 0, 1))
Coefficients:
         ar1      ma1  intercept
      0.8217  -0.5076    -0.2914
s.e.  0.1008   0.1530     0.2709
sigma^2 estimated as 1.019:  log likelihood = -142.99,  aic = 291.98
&/strong&&strong&
&/strong&&strong&&/strong&&p align=&left&&我们随机生成一些ARIMA模型的数据,并且假装我们不知道模型的参数,然后做估计练习。&/p&&p align=&left&& 生成&span style=&background-color: rgb(255, 204, 204);&&ARIMA数据的命令为arima.sim()&/span&. 自动回归的参数为0.9,移动平均的参数为-0.5,估计中要假设这两个参数的设置为未知,并进行估计。计算估计值的命令为&span style=&background-color: rgb(255, 204, 204);&&arima(&/span&)。其中,x指随机生成的数据。&/p&&p align=&left&&在下面的输出中,我们可以看到&span style=&background-color: rgb(255, 204, 204);&&AR参数的估计值为0.8465, ma的参数估计值为-0.5021&/span&,这两个估计值与前面生成数据时设置的参数非常相似。也就是说,估计得比较精确。&/p&
&/pre&&pre code_snippet_id=&485640& snippet_file_name=&blog__1175724& name=&code& class=&plain&&&strong style=&background-color: rgb(255, 204, 204);&&对估计的模型进行检测和评估,需要执行的命令为tsdiag()&/strong&&strong style=&background-color: rgb(255, 255, 255);&&tsdiag(x.fit, gof.lag=20) # you know the routine- ?tsdiag for details
x = arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=100) # simulate some data
(x.fit = arima(x, order = c(1, 0, 1)))   # fit the model and print the results
tsdiag(x.fit, gof.lag=20) # you know the routine- ?tsdiag for details&/strong&&strong style=&background-color: rgb(255, 255, 255);&&
&/strong&&span style=&background-color: rgb(255, 255, 255);&&&strong&分析:&/strong&&/span&&p align=&left&&&strong&&span style=&background-color: rgb(255, 255, 255);&&第一个图表代表估计模型误差的绘图。英文叫做Standardized Residuals, 上面有很多竖线在横向坐标的上下分布。如果这个估计的模型比较可信,竖线的长度是比较相似的。如果竖线的长度互相有很大出入或者根本就不同,估计模型的可信度就非常差。下面误差绘图中竖线的长度比较相似,都处在稳定范围之内,即估计的模型没产生不符合要求的误差分布。(运行的结果与书上的不大一样)&/span&&span style=&background-color: rgb(255, 204, 204);&&可根据散点的范围波动来判断&/span&&/strong&&/p&&p align=&left& style=&font-weight: background-color: rgb(255, 255, 255);&& &/p&&p align=&left& style=&font-weight: background-color: rgb(255, 255, 255);&&再介绍输出的第二张绘图,标题是ACF of Residuals。ACF指数据点相互之间的关系,当然在生成这个数据时,数据点之间互相独立,并不存在任何关系。所以在这张图上,只有位于0刻度上的竖线最高,其ACF值为1。这个0代表数据点与自己相比较, 即数据点永远和它自己有关系,这种关系数值为1。其他横向数轴上的刻度代表一个数据点于其他数据点之间的关系,这些刻度上竖线的长度几乎等于0,即这个数据点与其他数据点没明显关系。这张ACF图代表估计的模型没造成误差之间的任何关系。这是符合数据生成时每个数据都是独立的这个前提的。由此可见,这ACF图符合检测要求。&/p&&p align=&left& style=&font-weight: background-color: rgb(255, 255, 255);&& &/p&&p align=&left& style=&font-weight: background-color: rgb(255, 255, 255);&&下面来介绍第三张图,也就是Ljung-Box 指标。这个指标可对每一个时间序列的延迟进行显著性的评估。这张图的横坐标代表时间序列的延迟,纵坐标代表P-value,即显著性。如果P-value十分小,就说明在其相对应的延迟点上是显著的。我们就需要抛弃所假设的模型,并且结论在所假设的模型不可信。需要注意的是,他们使用假设的模型对一个时间序列进行估计,如果P-value是显著的话,我们使用的模型就不可信,需要尝试其他新模型。具体判定技巧是,P-value点的高度越高,我们的模型越可信。&/p&&strong style=&background-color: rgb(255, 255, 255);&&
&span style=&background-color: rgb(255, 204, 204);&&预测(predict)&/span&x = arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=100) # simulate some data
(x.fit = arima(x, order = c(1, 0, 1)))   # fit the model and print the results
x.fore = predict(x.fit, n.ahead=10)  
# plot the forecasts
U = x.fore$pred + 2*x.fore$se
L = x.fore$pred - 2*x.fore$se
minx=min(x,L)
maxx=max(x,U)
ts.plot(x,x.fore$pred,col=1:2, ylim=c(minx,maxx))
lines(U, col=&blue&, lty=&dashed&)
lines(L, col=&blue&, lty=&dashed&) 
&strong&分析:&/strong&&p align=&left&&&strong&预测所使用的语句是predict(). 在()中输入x.fit 和n.ahead =10 ,其中x.fit 代表前面估计出来的模型与参数,n.ahead =10 代表向前预测10个数据点。这条语句给我们生成的预测程序项目,被命名为x.fore. &/strong&&/p&&p align=&left&&&strong&在下面两条语句中,计算了两列的新数据,分别被称为U和L. U指最高界限,L指最低界限。 最高限和最低限告诉我们未来数据有很高可能性事发生在两个界限之间,超出或低于这两个界限的几率并不高。&/strong&&/p&&p align=&left&&&strong& 下一步使用Min 和 Max条语句,把数据中的最大值和最小值挑选出来,为接下来的统计绘图设置图标界限。再使用TS.PLOT()语句. 把原有的数据和向前预测10个数据点的数据绘在一张图上。最后使用lines()语句,把计算的最高界限和最低界限加到绘好的图上。这两条界限的颜色为蓝色,形式也被设置为虚线的形式。&/strong&&/p&
&strong&除去原有的数据外,又加上了三条曲线。最上面和最下面的曲线是虚线,是计算出的最高界限和最低界限。中间的曲线为实线,是对未来的10个数据点进行的预测。在计算最高、最低限时,我们使用了一个2的系数,也就是说,这两个界限之间的几率为95%。 在未来的10个数据点有95%的几率是出现在上限和下限之间的。而出现在两个界限之外的几率小于或等于5%。&/strong&
(4)漂移参数的估计方法
u =read.table(&http://www.stat.pitt.edu/stoffer/tsa2/data/globtemp2.dat&)
# read the data
gtemp = ts(u[,2], start=1880, freq=1)
# yearly temp in col 2
plot(gtemp)
arima(gtemp, order=c(1,1,1))
drift = 1:length(gtemp) #生成了一个名为drift的空数列,其长度与全球温度数列的长度相同
arima(gtemp, order=c(1,1,1), xreg=drift)
Coefficients:
生成了一个名为drift的空数列,其长度与全球温度数列的长度相同。而在第二个语句中,我们只要把drift数列作为xreg输入就可以了。从输出中,我们可以看到AR参数值为0.2695,MA参数值为-0.8180,是drift值为0.0061,即漂移产生的值. 上面提到过漂移值为+0.6,我们的估计值比提供的数值小100倍。其原因是因为我们使用的全球温度数据是经过被100整除的,所以估计的漂移值也相应地小了100倍。如果将gtemp数列输出后看一下,我们就可以明白所有温度都是被100整除过的。
(5)运用R语言研究非独立误差与线性回归
library(nlme) #load the package
trend = time(mort) #assumes mort and part are there from previous examples
fit.lm = lm(mort~trend + part)
acf(resid(fit.lm))
# check acf and pacf of the resids
pacf(resid(fit.lm))
# or use acf2(resid(fit.lm)) if you have acf2
&p align=&left&&pacf绘图是第五行命令运行的结果。在上图中,横向坐标代表的延迟,被叫做lag. 在lag1和lag2上出现了高于蓝色虚线的pacf,其他lag上超出蓝色虚线的距离并不明显可以忽略。这两个高于蓝色虚线的pacf说明这个误差数列符合AR2模型。&/p&
fit.gls = gls(mort~trend + part, correlation=corARMA(p=2), method=&ML&) # resids appear to be AR(2) ... now use gls() from nlme:使用gls()命令来对数据建立带有AR2误差的线性回归模型,生成结果记为fit.gls
summary(fit.gls)
(fit2.gls = arima(mort, order=c(2,0,0), xreg=cbind(trend, part)))
##xreg=cbind(trend,part)代表用时间和空气污染物作为进行预测的变量
Coefficients:
sigma^2 estimated as 28.99:
log likelihood = -1576.56,
aic = 3165.13
Box.test(resid(fit2.gls), 12, type=&Ljung&)
# and so on ... #Box.test计算n独立的null假设
(6)运用R语言研究估计波动的频率
x = arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=2^8) # some data
(u = polyroot(c(1,-1,.9)))
# x is AR(2) w/complex roots
Arg(u[1])/(2*pi)
# dominant frequency around .16:
par(mfcol=c(3,1))
plot.ts(x)
spec.pgram(x, spans=c(3,3), log=&no&)
# nonparametr also see spectrum()
#对已生成的随机波动数据的频率进行估计。注意,在这里不能够使用自动回归模型作为已知条件来进行频率的计算,而是要对频率进行估计。
spec.ar(x, log=&no&)
# parametric spectral estimate
17.&A Little Book of R for Time Series &Examples
(1)数据集
##########################################
#A Little Book of R for Time Series
##########################################
##########################################
#读取时间序列数据
##########################################
####################################################
#http://robjhyndman.com/tsdldata/misc/kings.dat
#包含着从威廉一世开始的英国国王的去世年龄数据
####################################################
####################################################
#读取时间序列数据
####################################################
kings &- scan(&http://robjhyndman.com/tsdldata/misc/kings.dat&, skip = 3)
####################################################
#将读入的数据存入到一个时间序列对象中
####################################################
kingstimeseries &- ts(kings)
kingstimeseries
####################################################
#ts可以用于指定frequency, start起始值,end结束值
####################################################
##########################################################
#一个样本数据集是从1946年1月到1959年12月的纽约每月
#出生人口数量(由牛顿最初收集)数据集可以从此链接下
#载(http://robjhyndman.com/tsdldata/data/nybirths.dat)
##########################################################
births &- scan(&http://robjhyndman.com/tsdldata/data/nybirths.dat&)
birthstimeseries &- ts(births, frequency = 12, start = c(1946, 1))
birthstimeseries
##########################################################
#http://robjhyndman.com/tsdldata/data/fancy.dat
#包含着一家位于昆士兰海滨度假圣地的纪念品商店
#从1987年1月到1987年12月的每月销售数据
##########################################################
souvenir &- scan(&http://robjhyndman.com/tsdldata/data/fancy.dat&)
souvenirtimeseries &- ts(souvenir, frequency = 12, start = c(1987, 1))
souvenirtimeseries
##########################################################
#plotting Time Series绘制时间序列图
##########################################################
##########################################################
#可使用R中的plot.ts()函数来画时间序列图
##########################################################
plot.ts(kingstimeseries)
##########################################################
#变化不大,可用相加模型来描述
##########################################################
##########################################################
#画出一个纽约每月出生人口数量的时间序列图
##########################################################
plot.ts(birthstimeseries)
##########################################################
#这个时间序列在一定月份存在的季节性变动:在每年的夏天都
#有一个出生峰值,在冬季的时候进入波谷。同样,这样的时间
#序列也可能是一个相加模型,随着时间推移,季节性波动时大
#致稳定的而不是依赖于时间序列水平,且对着时间的变化,随
#机波动看起来也是大致稳定的。
##########################################################
##########################################################
#画出澳大利亚昆士兰州海滨度假圣地的纪念品商店从1987年1月
#到1987年12月的每月销售数据。
##########################################################
plot.ts(souvenirtimeseries)
##########################################################
#在这个案例中,看上去似乎相加模型不适合描述这个时间序列,
#因为这个季节性波动和随机变动的大小是随着时间序列逐步上
#升的水平。因此,我们需要将时间序列进行转换,以便得到一
#个可以用相加模型描述的时间序列。例如,我们对原始数据取
#自然对数进行转换计算:
##########################################################
logsouvenirtimeseries &- log(souvenirtimeseries)
plot.ts(logsouvenirtimeseries)
##########################################################
#可以看到季节性波动和随机变动的大小在对数变换后的时间序列
#上,随着时间推移,季节性波动和随机波动的大小是大致恒定的,
#并且不依赖于时间序列水平。因此,这个对数变换后的时间序列
#也许可以用相加模型进行描述。
##########################################################
(2)分解时间序列
##########################################################
#分解时间序列
##########################################################
##########################################################
#分解一个时间序列意味着把它拆分成构成元件,一般序列包含
#一个趋势部分、一个不规则部分,如果是一个季节性时间序列,
#则还有一个季节性部分。
##########################################################
##########################################################
#分解非季节性数据
##########################################################
##########################################################
#一个非季节性时间序列包含一个趋势部分和一个不规则部分。
#分解时间序列即为试图把时间序列拆分成这些成分,也就是说,
#需要估计趋势的和不规则的这两个部分
#为了估计出一个非季节性时间序列的趋势部分,
#使之能够用相加模型进行描述,最常用的方法便是平滑法,
#比如计算时间序列的简单移动平均。
#在R的“TTR”包中的SMA()函数可以用
#简单的移动平均来平滑时间序列数据。
##########################################################
library(TTR)
##########################################################
#SMA用于计算序列中的不同的MA
##########################################################
##########################################################
#使用SMA()函数去平滑时间序列数据
#通过参数&n&指定来简单移动平均的
##########################################################
##########################################################
#kingstimeseries的plot呈现出非季节性
#且由于其随机变动在整个时间段内是大致不变的
#这个序列也可被描述为一个相加模型
#故可尝试使用简单移动平均平滑来估计趋势部分
#采用跨度不3的简单移动平均平滑时间序列数据
##########################################################
kingstimeseriesSMA3 &- SMA(kingstimeseries, n = 3)
plot.ts(kingstimeseriesSMA3)
##########################################################
#当我们使用跨度为3的简单移动平均平滑后,
#时间序列依然呈现出大量的随便波动。
#因此,为了更加准确地估计这个趋势部分,
#我们也许应该尝试下更大的跨度进行平滑。
#正确的跨度往往是在反复试错中获得的。
##########################################################
kingstimeseriesSMA8 &- SMA(kingstimeseries, n = 8)
plot.ts(kingstimeseriesSMA8)
##########################################################
#这个跨度为8的简单移动平均平滑数据的趋势部分看起来更加清晰了
#我们可以发现这个时间序列前20为国王去世年龄从最初的55周岁下
#降到38周岁,然后一直上升到第40届国王的73周岁。
##########################################################
##########################################################
#分解季节性数据
##########################################################
##########################################################
#一个季节性时间序列包含:
#一个趋势部分,一个季节性部分和一个不规则部分。
#分解时间序列就意味着要把时间序列分解称为这三个部分:
#也就是估计出这三个部分。
##########################################################
##########################################################
#对于可以使用相加模型进行描述的时间序列中的趋势部分和
#季节性部分,我们可以使用R中“decompose()”函数来估计。
#这个函数可以估计出时间序列中趋势的、季节性的和不规则
#的部分,而此时间序列须是可以用相加模型描述的。
##########################################################
##########################################################
#对于可以使用相加模型进行描述的时间序列中的趋势部分和
#季节性部分,我们可以使用R中“decompose()”函数来估计。
#这个函数可以估计出时间序列中趋势的、季节性的和不规则
#的部分,而此时间序列须是可以用相加模型描述的。
##########################################################
##########################################################
#“decompose()”这个函数返回的结果是一个列表对象,
#里面包含了估计出的季节性部分,趋势部分和不规则部分,
#他们分别对应的列表对象元素名为“seasonal”、“trend”、和“random”。
##########################################################
##########################################################
#为了估计时间序列的趋势的、季节性和不规则部分,输入代码:
##########################################################
birthstimeseriescomponents &- decompose(birthstimeseries)
birthstimeseriescomponents
##########################################################
#估计出的季节性、趋势的和不规则部分现在被存储在变量
#birthstimeseriescomponents$seasonal,
#birthstimeseriescomponents$trend 和
#birthstimeseriescomponents$random 中##########################################################
##########################################################
#可以使用“plot()”函数画出时间序列中估计的趋势的、
#季节性的和不规则的部分
##########################################################
plot(birthstimeseriescomponents)
##########################################################
#季节性因素调整
##########################################################
##########################################################
#如果这个季节性时间序列可以用相加模型来描述,
#你可以通过估计季节性部分修正时间序列,
#也可以从原始序列中去除掉估计得季节性部分。
#我们可以通过“decompose()”函数使用估计出的季节性部分进行计算。
##########################################################
##########################################################
#例如,对纽约每月出生人口数量进行季节性修正,
#我们可以用“decompose()”估计季节性部分,
#也可以把这个部分从原始时间序列中去除。
##########################################################
birthstimeseriescomponents &- decompose(birthstimeseries)
birthstimeseriesseasonallyadjusted &- birthstimeseries - birthstimeseriescomponents$seasonal
##########################################################
#可以使用“plot()”画出季节性修正时间序列,代码如下
##########################################################
plot(birthstimeseriesseasonallyadjusted)
##########################################################
#这个季节性修正后的时间序列现在仅包含趋势部分和不规则变动部分。
##########################################################
(3)使用指数平滑法进行预测
##########################################################
#使用指数平滑法进行预测
##########################################################
##########################################################
#指数平滑法可以用于时间序列数据的短期预测。
##########################################################
##########################################################
#简单指数平滑法
##########################################################
##########################################################
#如果你有一个可用相加模型描述的,
#并且处于恒定水平和没有季节性变动的时间序列
#指数平滑法可以用于时间序列数据的短期预测。
##########################################################
##########################################################
# http://robjhyndman.com/tsdldata/hurst/precip1.dat
#这个文件包含了伦敦从1813年到1912年全部的每年每英尺降雨量
#(初始数据来自Hipel and McLeod, 1994)
##########################################################
rain &- scan(&http://robjhyndman.com/tsdldata/hurst/precip1.dat&, skip = 1)
rainseries &- ts(rain, start = c(1813))
plot.ts(rainseries)
##########################################################
#从这个图可以看出整个曲线处于大致不变的水平
#(意思便是大约保持的25英尺左右)。
#其随机变动在整个时间序列的范围内也可以认为是大致不变的,
#所以这个序列也可以大致被描述成为一个相加模型。
#因此,我们可以使用简单指数平滑法对其进行预测。
##########################################################
##########################################################
#为了能够在R中使用简单指数平滑法进行预测,
#我们可以使用R中的“HoltWinters()”函数对预测模型进行修正。
#为了能够在指数平滑法中使用HotlWinters(),
#我们需要在HoltWinters()函数中设定参数beta=FALSE和gamma=FALSE
##########################################################
##########################################################
#HoltWinters()函数返回的是一个变量列表,包含了一些元素名。
##########################################################
rainseriesforecasts &- HoltWinters(rainseries, beta = FALSE, gamma = FALSE)
rainseriesforecasts
##########################################################
#Holt-Winters exponential smoothing without trend and without seasonal component.
HoltWinters(x = rainseries, beta = FALSE, gamma = FALSE)
#Smoothing parameters:
#beta : FALSE
#gamma: FALSE
#Coefficients:
#a 24.67819
##########################################################
##########################################################
#HoltWinters()的输出告诉我们alpha参数的估计值约为0.024。这个数字
#非常接近0,告诉我们预测是基于最近的和较远的一些观测值
#(尽管更多的权重在现在的观测值上)。
##########################################################
##########################################################
#我们将HoltWinters()函数的输出结果存储在
#“rainseriesforecasts”这个列表变量里。
#这个HoltWinters()产生的预测呗存储在一个元素名为“fitted”的列表变量里,
#我们可以通过以下代码获得这些值:
##########################################################
rainseriesforecasts$fitted
##########################################################
#可以再画出原始时间序列和预测的,代码如下:
##########################################################
plot(rainseriesforecasts)
##########################################################
#这个图用黑色画出了原始时间序列图,用红色画出了预测的线条。
#在这里,预测的时间序列比原始时间序列数据平滑非常多。
##########################################################
##########################################################
#作为预测准确度的一个度量,
#我们可以计算样本内预测误差的误差平方之和,
#即原始时间序列覆盖的时期内的预测误差。
#这个误差平方法将存储在一个元素名为
#“rainseriesforecasts”(我们称之为“SSE”)的列表变量里
##########################################################
rainseriesforecasts$SSE
##########################################################
#rainseriesforecasts$SSE
##########################################################
##########################################################
#用时间序列的第一个值作为这个水平的初始值在简单指数平滑法
#中常见的操作。例如,在伦敦降雨量这个时间序列里,第一个值
#为1813年的23.56(英尺)。你可以在HoltWinters()函数中使
#用“l.start”参数指定其味初始值。例如,我们将预测的初始值
#水平设定为23.56,代码如下:
##########################################################
HoltWinters(rainseries, beta = FALSE, gamma = FALSE, l.start = 23.56)
##########################################################
#可以使用R中的“forecast”包中的“forecast.HoltWinters()”
#函数进行更远时间点上的预测。使用Forecast.HoltWinters()函数,
#我们首先得安装R的“forecast”包
##########################################################
install.packages(&forecast&)
library(forecast)
##########################################################
#当我们使用forecast.HoltWinters()函数时,
#如它的第一个参数(input),
#你可以在已使用HoltWinters()函数调整后的预测模型中忽略它。
#例如,在下雨的时间序列中,
#使用HoltWinters()做成的预测模型存储在“rainseriesforecasts”变量中。
#你可以使用forecast.HoltWinters()中的参数”h”来制定你想要做多少时间点的预测
##########################################################
##########################################################
#要使用forecast.HoltWinters()做年(之后8年)的下雨量预测,我们输入:
##########################################################
rainseriesforecasts2 &- forecast.HoltWinters(rainseriesforecasts, h = 8)
rainseriesforecasts2
##########################################################
#Point Forecast
24.93 30.69 33.09470
24.33 30.24 33.09715
24.73 30.79 33.09960
24.13 30.34 33.10204
24.53 30.90 33.10449
24.94 30.45 33.10694
24.34 30.01 33.10938
24.74 30.56 33.11182
##########################################################
##########################################################
#orecast.HoltWinters()函数给出了一年的预测,
#一个80%的预测区间和一个95%的预测区间的两个预测。
##########################################################
plot.forecast(rainseriesforecasts2)
##########################################################
#这的蓝色线条是预测的的降雨量,
#深灰色阴影区域为80%的预测区间,
#淡灰色阴影区域为95%的预测区间。
##########################################################
##########################################################
#使用forecast.HoltWinters()返回的样本内预测误差将被存储
#在一个元素名为“residuals”的列表变量中。
#如果预测模型不可再被优化,连续预测中的预测误差是不相关的。
#换句话说,如果连续预测中的误差是相关的,
#很有可能是简单指数平滑预测可以被另一种预测技术优化。
##########################################################
##########################################################
#为了验证是否如此,我们获取样本误差中1-20阶的相关图。
#我们可以通过R里的“acf()”函数计算预测误差的相关图。
#为了指定我们想要看到的最大阶数,可以使用acf()中的“lag.max”参数。
##########################################################
acf(rainseriesforecasts2$residuals, lag.max = 20)
##########################################################
#可以从样本相关图中看出自相关系数在3阶的时候触及了置信界限。
#为了验证在滞后1-20阶(lags 1-20)时候的非0相关是否显著,
#我们可以使用Ljung-Box检验。这可以通过R中的“Box.test()”函数实现。
#最大阶数我们可以通过Box.test()函数中的“lag”参数来指定。
##########################################################
Box.test(rainseriesforecasts2$residuals, lag = 20, type = 'Ljung-Box')
##########################################################
#Box-Ljung test
rainseriesforecasts2$residuals
#X-squared = 17.4008, df = 20, p-value = 0.6268
##########################################################
##########################################################
#这里Ljung-Box检验统计量为17.4,并且P值是0.6,所以这是不足
#以证明样本内预测误差在1-20阶是非零自相关的。
#为了确定预测模型不可继续优化,
##我们需要一个好的方法来检验预测误差是正态分布,
#并且均值为零,方差不变。为了检验预测误差是方差不变的,
#我们可以画一个样本内预测误差图:
##########################################################
plot.ts(rainseriesforecasts2$residuals)
##########################################################
##为了检验预测误差是均值为零的正太分布,
#我们可以画出预测误差的直方图,
#并覆盖上均值为零、标准方差的正态分布的曲线图到预测误差上
##########################################################
plotForecastErrors &- function(forecasterrors){
# make a red histogram of the forecast errors:
mybinsize &- IQR(forecasterrors)/4
mysd &- sd(forecasterrors)
mymin &- min(forecasterrors) + mysd*5
mymax &- max(forecasterrors) + mysd*3
mybins &- seq(mymin, mymax, mybinsize)
hist(forecasterrors, col=&red&, freq=FALSE, breaks=mybins) # freq=FALSE ensures the area under the histogram = 1
# generate normally distributed data with mean 0 and 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 of the histogram of forecast errors:
points(myhist$mids, myhist$density, type=&l&, col=&blue&, lwd=2)
plotForecastErrors(rainseriesforecasts2$residuals)
##########################################################
#图展现出预测误差大致集中分布在零附近,
#或多或少的接近正太分布,
#尽管图形看起来是一个偏向右侧的正态分布。
#然后,右偏是相对较小的
#,我们可以可以认为预测误差是服从均值为零的正态分布。
##########################################################
(4)霍尔特指数平滑法
##########################################################
#霍尔特指数平滑法
##########################################################
##########################################################
#霍尔特指数平滑法
#时间序列可以被描述为一个增长或降低趋势的、没有季节性的相加模型
#可以使用霍尔特指数平滑法对其进行短期预测。
##########################################################
##########################################################
#Holt指数平滑法估计当前时间点的水平和斜率。其平滑化是
#由两个参数控制的,alpha,用于估计当前时间点的水平,
#beta,用于估计当前时间点趋势部分的斜率。
#正如简单指数平滑法一样,alpha和beta参数都介于0到1之间,
#并且当参数越接近0,大多数近期的观测则将占据预测更小的权重。
##########################################################
##########################################################
#数据来源:
#可能可以用相加模型描述的有趋势的、无季节性的时间序列案例就是这
#1866年到1911年每年女人们裙子的直径。这个数据可以从该文件获得
# http://robjhyndman.com/tsdldata/roberts/skirts.dat
#初始数据来自Hipel and McLeod, 1994
##########################################################
skirts &- scan(&http://robjhyndman.com/tsdldata/roberts/skirts.dat&, skip = 5)
skirtsseries &- ts(skirts, start = c(1866))
plot.ts(skirtsseries)
##########################################################
#可以从此图看出裙子直径长度从1866年的600增加到,
#并且在此之后有下降到1911年的520。
#为了进行预测,我们使用R中的HoltWinters()函数对预测模型进行调整。
#为了使用HoltWinters()进行Holt指数平滑法,
#我们需要设定其参数gamma=FALSE(gamma参数常常用于Holt-Winters指数平滑法)
##########################################################
skirtsseriesforexasts &- HoltWinters(skirtsseries, gamma = FALSE)
##########################################################
#Holt-Winters exponential smoothing with trend and without seasonal component.
# HoltWinters(x = skirtsseries, gamma = FALSE)
#Smoothing parameters:
alpha: 0.8383481
#gamma: FALSE
#Coefficients:
#a 529.308585
##########################################################
skirtsseriesforexasts$SSE
##########################################################
# skirtsseriesforexasts$SSE
#[1] 16954.18
##########################################################
##########################################################
#这里的alpha预测值为0.84,beta预测值为1.00。
#这都是非常高的值,告诉我们无论是水平上,还是趋势的斜率,
#当前值大部分都基于时间序列上最近的观测值。这样的直观感觉很好,
#因为其时间序列上的水平和斜率在整个时间段发生了巨大的变化。
#预测样本内误差的误差平方和是16954。
##########################################################
plot(skirtsseriesforexasts)
##########################################################
#从该图我们可以看到样本内预测非常接近观测值,
#尽管他们对观测值来说有一点点延迟。
#可以通过HoltWinters()函数中的“l.start”和“b.start”参数
#去指定水平和趋势的斜率的初始值。常见的设定水平初始值是
#让其等于时间序列的第一个值(在裙子数据中是608),而斜率的初始值则是其第
#二值减去第一个值(在裙子数据中是9)。例如,为了使用Holt指数平
#滑法找到一个在裙边直径数据中合适的预测模型,我们设定其水平初始值为608,趋势部分的斜率初始值为9
##########################################################
HoltWinters(skirtsseries, gamma = FALSE, l.start = 608, b.start = 9)
##########################################################
#Holt-Winters exponential smoothing with trend and without seasonal component.
HoltWinters(x = skirtsseries, gamma = FALSE, l.start = 608, b.start = 9)
#Smoothing parameters:
# alpha: 0.8346775
#gamma: FALSE
#Coefficients:
#a 529.278637
##########################################################
##########################################################
#可以使用“forecast”包中的forecast.HoltWinters()函数预测未来时间而无需覆盖原始序列
#现在有的1866年到1911年的裙边直径时间序列数据,
#因此我们可以预测1912年到1930年(19个点或者更多),并且画出他们
##########################################################
skirtsseriesfoecasts2 &- forecast.HoltWinters(skirtsseriesforexasts, h = 19)
plot.forecast(skirtsseriesfoecasts2)
##########################################################
#预测的部分使用蓝色的线条标识出来了,
#浅蓝色阴影区域为80%预测区间,灰色阴影区间为95%的预测区间。
##########################################################
##########################################################
#简单指数平滑法一样,我们瞧瞧样本内预测误差是否在延迟1-20阶时是非零自相关的,
#以此来检验模型是否还可以被优化
#可以创建一个相关图,进行Ljung-Box检验,代码如下:
##########################################################
acf(skirtsseriesfoecasts2$residuals, lag.max = 20)
Box.test(skirtsseriesfoecasts2$residuals, lag = 20, type = 'Ljung-Box')
##########################################################
#Box-Ljung test
skirtsseriesfoecasts2$residuals
#X-squared = 19.7312, df = 20, p-value = 0.4749
##########################################################
##########################################################
#这个相关图呈现出样本内预测误差的样本自相关系数在#
#滞后5阶的时候超过了置信边界。不管怎样,我们可以界
#定在前20滞后期中有1/20的自相关值超出95%的显著边界
#是偶然的,当我们进行Ljung-Box检验时,P值为0.47,
#意味着我们是不足以证明样本内预测误差在滞后1-20阶的
#时候是非零自相关的。
##########################################################
##########################################################
#和简单指数平滑法一样,我们应该检查整个序列中的预测误差是
#否是方差不变、服从零均值正态分布的。我们可以画出一个时间
#段预测误差图,和一个附上正太曲线的预测误差分布的直方图:
##########################################################
plot.ts(skirtsseriesfoecasts2$residuals)
plotForecastErrors(skirtsseriesfoecasts2$residuals)
##########################################################
#预测误差的时间曲线图告诉我们预测误差在整个时间段内是大致方差不变的
#这个预测误差的直方图告诉我们预测误差似乎是零均值、方差不变的正态分布。
#Ljung-Box检验告诉我们这是不足以证明预测误差是自相关的,
#而其预测误差的时间曲线图和直方图表示出似乎预测误差是服从零均值、方差不变的正态分布的。
#因此,我们可以总结这Holt指数平滑法为裙边直径提供了一个合适的预测,并且是不可再优化的
#另外,这也意味着基于80%预测区间和95%预测区间的假设是非常合理的。
##########################################################
(5)Holt-Winters指数平滑法
##########################################################
#Holt-Winters指数平滑法
##########################################################
##########################################################
#有一个增长或降低趋势
#并存在季节性
#可被描述成为相加模型的时间序列,
#可以使用霍尔特-温特指数平滑法对其进行短期预测。
##########################################################
##########################################################
#Holt-Winters指数平滑法估计当前时间点的水平,斜率和季节性部分。
#平滑化依靠三个参数来控制:alpha,beta和gamma,
#分别对应当前时间点上的水平,趋势部分的斜率和季节性部分。
#参数alpha,beta和gamma的取值都在0和1之间,
#并且当其取值越接近0意味着对未来的预测值而言最近的观测值占据相对较小的权重。
##########################################################
##########################################################
#用相加模型描述的并附有趋势性和季节性的时间序列案例,
#便是澳大利亚昆士兰州的海滨纪念品商店的月度销售日志
#可以使用HoltWinters()函数对预测模型进行修正
##########################################################
logsouvenirtimeseries &- log(souvenirtimeseries)
souvenirtimeseriesforecasts &- HoltWinters(logsouvenirtimeseries)
souvenirtimeseriesforecasts
##########################################################
#Holt-Winters exponential smoothing with trend and additive seasonal component.
# HoltWinters(x = logsouvenirtimeseries)
#Smoothing parameters:
alpha: 0.413418
#gamma: 0.9561275
#Coefficients:
##########################################################
##########################################################
#这里alpha,beta和gamma的估计值分别是0.41,0.00和0.96。
#alpha(0.41)是相对较低的,意味着在当前时间点估计得
#水平是基于最近观测和历史观测值。beta的估计值是0.00,
#表明估计出来的趋势部分的斜率在整个时间序列上是不变的,
#并且应该是等于其初始值。这是很直观的感觉,水平改变非常多,
#但是趋势部分的斜率b却仍然是大致相同的。与此相反的,
#gamma的值(0.96)则很高,表明当前时间点的季节性部分的估计仅仅基于最近的观测值。
##########################################################
plot(souvenirtimeseriesforecasts)
##########################################################
#可以从图中看出Holt-Winters指数平滑法是非常成功得预测了季节峰值,
#其峰值大约发生在每年的11月份。
##########################################################
##########################################################
#为了预测非原始时间序列的未来一段时间,
#我们使用“forecast”包中的“forecast.HoltWinters()”函数
#例如,纪念品销售的原始数据是1987年1月到1993年12月。
#如果我们想预测1994年1月到1998年12月(48月或者更多),并且画出预测,代码如下:
##########################################################
souvenirtimeseriesforecasts2 &- forecast.HoltWinters(souvenirtimeseriesforecasts, h = 48)
plot.forecast(souvenirtimeseriesforecasts2)
##########################################################
#蓝色线条显示出来的是预测,灰蓝和灰色阴影分别是80%和95%的预测区间。
##########################################################
&/pre&&p&&/p&&p&&/p&&pre code_snippet_id=&485640& snippet_file_name=&blog__5292070& name=&code& class=&plain&&##########################################################
#可以通过画相关图和进行Ljung-Box检验来检查样本内预测误差
#在延迟1-20阶时否是非零自相关的,并以此确定预测模型是否可以再被优化。
##########################################################
acf(souvenirtimeseriesforecasts2$residuals, lag.max = 20)
Box.test(souvenirtimeseriesforecasts2$residuals, lag = 20, type = 'Ljung-Box')
##########################################################
Box-Ljung test
souvenirtimeseriesforecasts2$residuals
#X-squared = 17.5304, df = 20, p-value = 0.6183
##########################################################
##########################################################
# 这个样本内预测误差的相关图并没有在延迟1-20阶内自相关系数超过置信界限的。
#而且,Ljung-Box检验的P值是0.6,意味着是不足以证明延迟1-20阶是非零自相关的。
#可以在整个时间段内检验预测误差是否是方差不变,并且服从零均值正态分布的。
#方法是画出预测误差的时间曲线图和直方图(并覆盖上正太曲线):
##########################################################
plot.ts(souvenirtimeseriesforecasts2$residuals)
plotForecastErrors(souvenirtimeseriesforecasts2$residuals)
##########################################################
#似乎告诉我们预测误差在整个时间段是方差不变的。
#从预测误差的直方图,似乎其预测误差是服从均值为零的正态分布的。
#这是不足以证明预测误差 在延迟 1-20 阶是自相关的,
#并且预测误差 自相关的,整个时间段呈现出服从零均值、方差 零均值、方差
#不变的正态分布。这 不变的正态分布。这 不变的正态分布。这 不变的正态分布。
#这暗示着 暗示着 Holt Holt-WintersWintersWinters WintersWintersWinters指数平滑法为
#纪念品商店的销售数据提供了一个合适预测模纪念品商店的销售数据提供了一个合适预测模型
#区间的假设 也是 合理 的。
##########################################################
(6)时间序列的差分
##########################################################
#对于平稳性正式的检验
#称作”单位根测试“
#可在fUnitRoots包中得到
##########################################################
##########################################################
#时间序列的差分
##########################################################
##########################################################
#ARIMAARIMAARIMAARIMAARIMA模型为平稳时间序列定义的。
#因此,如果你从一个非开始首先就需要做差分直 模型为平稳时间序列定义的。
#若必须对时间序列做d阶养分才能得到一个平稳序列,可使用ARIMA(p,d,q)模型,共中,d是差分的阶数
#可使用diff()函数作时间序列的差分
##########################################################
##########################################################
#我们可以通过键入下面的代码来得到时间序列(数据存于“skirtsseries“)
#的一阶养分,并画出差分序列的图
##########################################################
skirtsseriesdiff1 &- diff(skirtsseries, differences = 1)
plot.ts(skirtsseriesdiff1)
##########################################################
#一阶差分时间序列结果(上图) 一阶差分时间序结果均值看起来并不平稳
#因此,我们需要再次做差分来看一下是否能得到个平稳时间序列
##########################################################
skirtsseriesdiff2 &- diff(skirtsseries, differences = 2)
plot.ts(skirtsseriesdiff2)
##################################################}

我要回帖

更多关于 时间序列找趋势项和季节性的R代码 的文章

更多推荐

版权声明:文章内容来源于网络,版权归原作者所有,如有侵权请点击这里与我们联系,我们将及时删除。

点击添加站长微信