R语言时间序列中文教程可编辑.docx

R语言时间序列中文教程可编辑.docx

《R语言时间序列中文教程可编辑.docx》由会员分享，可在线阅读，更多相关《R语言时间序列中文教程可编辑.docx（29页珍藏版）》请在冰点文库上搜索。

R语言时间序列中文教程可编辑

R语言时间序列中文教程（可编辑）

R语言时间序列中文教程

R语言时间序列中文教程李智在乔治梅森2012特别声明:

R语言是免费语言，其代码不带任何质量保证，使用R语言所产生的后果由使用者负全责。

前言R语言是一种数据分析语言，它是科学的免费的数据分析语言，是凝聚了众多研究人员心血的成熟的使用范围广泛全面的语言，也是学习者能较快受益的语言。

在R语言出现之前，数据分析的编程语言是SAS。

当时SAS的功能比较有限。

在贝尔实验室里，有一群科学家讨论提到，他们研究过程中需要用到数据分析软件。

SAS的局限也限制了他们的研究。

于是他们想，我们贝尔实验室的研究历史要比SAS长好几倍，技术力量也比SAS强好几倍，且贝尔实验室里并不缺乏训练有素的专业编程人员，那么，我们贝尔实验室为什么不自己编写数据分析语言，来满足我们应用中所需要的特殊要求呢,于是，贝尔实验室研究出了S-PLUS语言。

后来，新西兰奥克兰大学的两位教授非常青睐S-PLUS的广泛性能。

他们决定重新编写与S-PLUS相似的语言，并且使之免费，提供给全世界所有相关研究人员使用。

于是，在这两位教授努力下，一种叫做R的语言在奥克兰大学诞生了。

R基本上是S-PLUS的翻版，但R是免费的语言，所有编程研究人员都可以对R语言做出贡献，且他们已经将大量研究成果写成了R命令或脚本，因而R语言的功能比较强大，比较全面。

研究人员可免费使用R语言，可通过阅读R语言脚本源代码，学习其他人的研究成果。

笔者曾有幸在奥克兰大学受过几年熏陶，曾经向一位统计系的老师提请教过一个数据模拟方面的问题。

那位老师只用一行R语句就解答了。

R语言的强大功能非常令人惊讶。

为

了进一步推广R语言，为了方便更多研究人员学习使用R语言，我们收集了R语言时间序列分析实例，以供大家了解和学习使用。

当然，这是非常简单的模仿练习，具体操作是，用复制粘贴把本材料中R代码放入R的编程环境;材料中蓝色背景的内容是相关代码和相应输出结果。

经过反复模仿，学习者便能熟悉和学会。

需要提醒学习者的是:

建议学习者安装了R语言编程，再继续阅读本材料;执行R命令时，请删除命令的中文注解，没使用过在命令中加入中文;如果学习者是初次接触R或者Splus，建议先阅读，如果学习者比较熟悉R语言，还可以阅读优秀时间序列读物EcomometricsinR，也可以上QuickR网站。

目录R语言时间序列分析1前言1目录21.运用R语言研究JJ数据32.运用R语言研究空气污染193.运用R语言研究自动回归移动平均集成模型224.运用R语言研究全球变暖理论315.运用R语言研究非独立误差与线性回归326.运用R语言研究估计波动的频率367.运用R语言研究厄尔尼诺频率398.运用R语言研究太阳黑子周期频率401.运用R语言研究JJ数据学习R言语时间序列分析程序操作，需要从最基础、最简单的学起，例如在命令窗口中，输入并执行2+2等于4的R语言命令。

2+2[1]4执行完2+2等于4的R语言命令后，我们可以开始时间序列，即学着把玩johnson&Johnson数据。

下载jj.dat或执行下面语句。

这个数据已被人上传到因特网中。

R所需要做的只是将网址进行扫描就可以将数据读取进入R的编程环境中。

下面有3种不同读取数据的方法:

jjscan"

#readthedata读取数据jj-scan"#readthedataanotherway

第二种方法读取数据scan"-jj#andanother第三种方法读取数据使用R语言的人，有的喜欢使用[-]，有的喜欢使用[-]，大多数医疗系统的工作者喜欢用[]，正因为如此才用了上面种不同读取数据的方法。

读取数据后，

键入并执行jj，数据在窗口便会有如下显示:

jj[1]0.710.630.850.44[5]0.610.690.92

0.55..........[77]14.04

12.9614.859.99[81]16.2014.6716.0211.61jj中有84个数据被称作

对象。

下面命令可以显示所有对象。

objects如果使用matlab，你会认为jj

是一个84行1列的向量，但实际上不是这样。

jj有次序，有长度，但没维度，

R称这些对象为向量，要小心区别。

在R里，矩阵有维度，但向量没维度。

这都

是程序语言的一些概念。

jj[1]#thefirstelement列中第一个数据[1]

0.71jj[84]#thelastelement列中最后一个数据[1]11.61jj[1:

4]#thefirst4elements列中第一至第四个数据[1]0.710.630.850.44jj[-1:

80]#everythingEXCEPTthefirst80elements列中除第80个以外的所

有数据[1]16.2014.6716.0211.61lengthjj#thenumberofelements

有多少个数据[1]84dimjj#butnodimensions...但没维度NULLnrowjj#...norows没行NULLncoljj#...andnocolumns没列NULL#

如果你要把jj转变为一个向量，执行如下操作后，维度为84行1列。

jj

as.matrixjjdimjj[1]841然后把jj转变为一个时间序列对象。

jjts

jj,start1960,frequency4#ts命令这个数据是从1960年开始，个个季度

的收入，frequency4指四个季度。

R语言的优势在于可用一条命令做很多事，

即可以把前面的命令放在一起打包执行。

其操作如下:

jjtsscan"",start1960,frequency4在上面命令里，scan可以被read.table替代。

用read.table

读取数据可生成matrix对象，还可以给每列起名字。

下面学习一下read.table,dataframes,和时间序列对象。

输入命令后，窗口会有如下显示:

jjts

read.table",start1960,frequency4helpread.tablehelptshelpdata.frame需要注意的是，Scan和read.table不一样。

Scan生成的是有维

度的向量，read.table生成的则是带有维度的数据架构。

读取jj数据的最后

要领。

如果这个数据是从1960年第三个季度开始的，所需输入命令则为ts

x,startc1960,3,frequency4;如果是一个每月每月的数据，例如数据是从

1984年6月开始的，需要输入的命令则为tsx,startc1984,6,frequency12。

输入命令后，转变后的时间序列对象为:

jjQtr1Qtr2Qtr3Qtr419600.71

0.630.850.4419610.610.690.92

0.55..........197914.0412.9614.859.99198016.2014.6716.0211.61注意到区别了吗,时间信

息，也就是4个不同的季度的数据被加载到里面了。

进行时间数据分析后，窗

口会有如下显示:

timejjQtr1Qtr2Qtr3Qtr419601960.001960.251960.501960.7519611961.001961.251961.50

1961.75............19791979.001979.251979.501979.7519801980.001980.251980.501980.75接下来输入如下组合命令。

jjtsscan",start1960,frequency4然

后进行对数据绘图:

plotjj,ylab"EarningsperShare",main"J&J"输

入以上命令后，可以看到如下结果:

再输入下面的命令，看看区别。

plotjj,type

"o",col"blue",lty"dashed"plotdifflogjj,main"loggedanddiffed"下面利用操作两个相关命令，显示区别。

x-5:

5#

sequenceofintegersfrom-5to5y5*cosx#guessparmfrowc3,2#multifiguresetup:

3rows,2cols#---plot:

plot

x,main"plotx"plotx,y,main"plotx,y"#---plot.ts:

plot.tsx,main"plot.tsx"plot.tsx,y,main"plot.tsx,y"#---ts.plot:

ts.plotx,main"ts.plotx"ts.plottsx,tsy,col1:

2,main"ts.plotx,y"#note-xandyaretsobjects#---thehelpfiles[?

andhelparethesame]:

?

plot.tshelpts.plot?

par#mightaswellskimthegraphicalparametershelpfilewhileyou'rehere从窗口中的显示可以看出，如果数

据是时间序列对象，使用plot命令就足够了;如果数据是平常序列，使用

plot.ts也可以做时间绘图。

不过，把jj数据放在一张图上，数据会随着时

间的变化上上下下跳动，能从整体上反应上升或者下降的趋势。

上文中用红色光

滑的曲线代表上升的趋势，简单明了。

这需要将过滤和光滑的技巧使用在jj数

据上。

在这里，我们用对称的移动平均值来达到过滤和光滑的目的。

下面使用公

式:

fjjt?

jjt-2+?

jjt-1+?

jjt+?

jjt+1+?

jjt+2除

此之外，lowess的过滤平滑技巧（蓝色曲线）也要使用在jj数据中。

具体操作

如下图:

kc.5,1,1,1,.5#kisthevectorofweights用

于对称移动平均的系数kk/sumk[1]0.1250.2500.2500.2500.125fjjfilterjj,sides2,k#?

filterforhelp[butyouknewthatalready]使用对称移动平均plotjjlinesfjj,col"red"#addsalinetotheexistingplot称移动平均的绘图lineslowessjj,col"blue",lty

"dashed"#lowess的绘图操作后，窗口会显示下面结果:

看完jj数据，我们

就需要开始具体分析。

第一步，我们把所有jj数据都取log值。

第二步，我们

把log值做差，即使用log值数列中第二值减去第一值，第三值减去第二值，第

四值减去第三值等等。

如果做差处理前数列里有n个数值，处理后的结果中将有

n-1个数值。

dljjdifflogjj#differencetheloggeddata做log和差的处理plotdljj#plotitifyouhaven'talready对结果制图shapiro.testdljj#testfornormality测

试结果的正态分布的性质Shapiro-Wilknormalitytestdata:

dljjW0.9725,p-value0.07211处理完毕以上两步，我们接下来就要将柱形分布图

和图放在一起。

这两个图的本质仍旧在于测试数据正态分布的性质。

数据正态分

布的性质是整个统计学构架坚实的基础，如果这个性质的存在比较可信、通过了

所有测试，统计分析中得出的结论就比较可信、就通得过所有测试。

当然如果这

个性质在数据中不存在，我们需要用其它的技巧来处理。

详细的，参看R语言样

品比较应用的实例。

以上操作，在窗口中有如下显示:

parmfrowc2,1#setupthegraphics设置为两图的输出histdljj,probTRUE,12#histogram柱形分布图linesdensitydljj#smoothit-?

densityfordetails柱形分布图的曲线normdljj#normalQ-Qplot图

linedljj#addaline在图上加直线经过测试数据后，

窗口会有如下显示:

在实践操作中，时间序列数据存在着前后关系。

例如，今天

股票的价格很有可能决定明天股票的价格。

明天的温度取决于今天的气温。

做天

气预报的具体操作方法，是使用已经存在的天气历史记录，比如说今天的气温，

昨天的气温，前天的气温等等，来预测明天的气温。

当然，在进行预测之前，我

们一定要看清时间序列数据中的前后关系结构，清楚哪一个特定的历史数据可以

精确预测未来的数据。

在这里，我们使用被log和求差后的dljj数据，来介绍

分析时间序列数据前后关系结构的具体技巧。

在预测的实际应用中，我们总希

望用历史数据来预测即将要产生的数据。

事实上，已产生的数据相对于即将产生

的数据，中间存在着一定的延迟，也就是lag.比方说在某天凌晨12点开始记

录室内温度，每小时记一次，一共连续记录了10个小时。

凌晨12点的数据和凌

晨3点的数据之间就存在着延迟。

12点的数据比3点的早了3个小时，可记作

lag3.3点的数据比12点的晚了3个小时，可记作lead3.我们下面来介绍关

联性。

例如，冷饮的销量与天气温度存在关联性。

温度越高冷饮销量就越高，温

度越低冷饮销量也越低。

这种关联性称为正面关联性。

又如，人的体重和跑步速

度也存在关联性。

不过，人的体重越重，跑步速度就会越慢，体重越低，相对来

讲，速度就会越快。

这种关联性称为负面关联性。

下面我们回到预测应用上。

如

果现在收集的数据与将来的数据之间存在着正面或者是负面的关联性，我们就可

以用现在收集的数据来预测未来的数据。

因此找到现在收集到的数据与未来数据

之间的关联性是最关键的。

找到这种关联性的具体技巧被称作延迟图表，也就是

lag.plot.我们可以在电脑里输入如下命令:

lag.plotdljj,9,do.linesFALSE

#whythedo.linesFALSE?

...tryleavingitout上面语句显示了lag.plot

用法，输入的数据是被log和作差后的jj.dat数据。

其中9表示我们要考虑从

1到9这9个不同的延迟。

值得注意的是，在下面图表中显示出延迟4和8显示

出了正面关系。

其他几个延迟存在着负面关系。

我们可以利用这4和8的正面关

系来进行预测，即用现有数据计算接下来的第4个或者是第8个数据。

下面我

们来看ACF和PACF图表。

确定我们已经观察到的正面和负面关系。

parmfrowc

2,1#Thepowerofaccurateobservationiscommonlycalledcynicism#bythosewhohavenotgotit.-GeorgeBernardShawacfdljj,20#ACFtolag20-nographshown...keepreadingpacfdljj,20#PACFtolag20-nographshown...keepreading#!

!

NOTE!

!

acf2onthelinebelow

isNOTavailableinR...detailsfollowthegraphbelowacf2dljj#

thisiswhatyou'llseebelow在上图中，ACF和PACF横坐标的标记是1,2,3,4,5.但因为数据是季度性的，每年有4个季度所以1,2,3,4,5的标记代表的4,8,12,16,20的延迟。

当然，如果我们不喜欢上面横坐标的标记，也可以将dljj更改为tsdljj,freq1;也就是说acftsdljj,freq1,20。

因为在上面ACF图表中横坐标1代表的是延迟4，横坐标2代表的是延迟6，其相应的竖线代表的就是延迟4和8的正面关系。

接下来，下面我们介绍结构拆析。

在前面R代码中，我们曾将所有jj数据进行了lag变型。

在变型后的数据中，存在着上升趋势，季节的影响和每一时间点产生的随机的误差。

根据这一数据图，我们能够把趋势、季节和误差从变型后的jj数据中拆析出来，分别研究，或者分别进行绘图，以便于单独检查。

下面等式代表将要使用的数学模型:

Logjj趋势+季节+误差logjjtrend+season+error结构拆析的R命令是stl,下面语句中stl命令中输入的是lag变型后的jj数据。

其中的“per”输入指的是使用季节循环来进行拆析。

stl语句在这里生成了一个叫dog的R物件，然后Plot语句将dog物件进行绘图。

具体操作如下图”plotdog-stllogjj,

"per"窗口会出现下面所示:

上图中有四行R的绘图，其中第一行代表原来的log（jj）的数据。

此数据可以看到总体的上升趋势还存在着一定季节循环性的变化。

第二行绘图代表拆析后季节循环的作用。

第三行绘图代表拆析后将季节循环作用清除剩余的上升趋势，此数据清楚地看到那种循环性变化已经不存在，剩余的只是趋势。

第四行绘图代表将季节循环作用和总体的趋势从数据中清除后所剩余的随机产生的误差。

如果我们需要对数据的误差进行一些常规检测，例如进行正态分布检测，绘制图，还有绘制柱形图。

我们所需要的具体误差数据被存

在叫做dog$time.series[,3]的数列里。

即叫dog的物件中有个叫time.series的数据矩阵，误差就被存储在这个数据矩阵的第三列里。

$指调取dog物件中的time.series数据矩阵。

[,3]指数据矩阵中第三列。

如果要对这一数列的误差值进行ACF的分析，只需要执行命令acfdog$time.series[,3]。

再接下来，我们对log变型后的jj数据进行线性回归模型分析。

与上面结构拆析不同的是，我们在这里使用四个季度来量化季节循环对数据的影响。

一年中有四个季度，也是我们所使用数据所代表的。

这个jj数据是某一家公司的季度收入数据，从上面绘图中我们就可以看到，每一年第三季度就会出现一个收入高峰，随之而来第四季度收入就会跌入低谷。

然后在一季度和二季度收入又会逐渐上升。

这也就是说，每一季度对这家公司收入的影响都是不一样的。

具体考虑到这种季度之间的不同，我们可使用如下数学模型:

logjjβ*time+α1*Q1+α2*Q2+α3*Q3+α4*Q4+ε这个数学模型的意思是:

logjj趋势*时间+一季度的影响*一季度+二季度的影响*二季度+三季度的影响*三季度+四季度的影响*四季度+误差上面的模型β代表的就是总体上升趋势，α1α2α3α4代表的是四个季度的影响。

有一个非常有趣的问题，上面模型是把所有四个季度的趋势都加在了一起，其结果却是某单一季度的收入。

四个季度的和如何能够与一个季度相等问题就出在Q1Q2Q3Q4上。

因为Q被我们称作指示性函数。

函数的意思就是数据进，数据出，也就是说把一个数据输入到一个函数中，那个函数就会输出一个结果。

以上面的Q1函数为例，Q1只能输出两种结果，1和0.Q1所需要的输入是四种1，2，3，4，代表四个季度。

把1输入到Q1函数中时，Q1函数输出的结果为1，当把2,3,4输入Q1函数时，Q1函数输出的结果为0.与Q1函数类似，Q2函数的输入也是1,2,3,4，但只有输入为2时，Q2函数的输出才为

1，当输入为1,3,4时，Q2函数的输出为0.Q3函数输入为1,2,3,4，只有当输

入为3时，输出为1，输入其他数据时，输出为0.Q4函数的输入为1,2,3,4，只

有当输入为4时，输出为1，其他数据时，输出为0.我们再回到上面的模型，当

一个数据是从第一季度中记录下来的，Q1给出数值1，Q2给出数值0，Q3给出

数值0，Q4给出的数值0。

因为这时Q2，Q3，Q4都是0，二季度，三季度，四季

度的影响被0相乘后也变成了0.所以在第一季度Q1为1，而其他的为0.我们

就只考虑了一季度的影响，其他季度的影响不存在。

同理，当季度为二、三、四

时也有类似结果。

下面是建立这个线性模型的R语句，只有头三行是用来生成线

性模型的，第四条语句summary用来输出模型参数数值。

具体操作以及显示如

下:

Qfactorrep1:

4,21#makeQuarterfactors[that'srepeat1,2,3,4,21times]trendtimejj-1970#notnecessaryto"center"time,buttheresultslooknicerreglmlogjj~0+trend+Q,na.actionNULL#runtheregressionwithoutanintercept#--thena.actionstatementistoretaintimeseriesattributessummaryregCall:

lmformulalogjj~0+trend+Q,na.actionNULLResiduals:

Min1QMedian3Q-0.29318-0.09062-0.011800.084600.27644Coefficients:

EstimateStd.ErrortvaluePr|t|trend0.1671720.00225974.002e-16***Q11.0527930.02735938.482e-16***Q21.0809160.02736539.502e-16***Q31.1510240.02738342.032e-16***Q40.8822660.02741232.192e-16***---Signif.codes:

0'***'0.001'**'0.01'*'0.05'.'0.1''1Residualstandarderror:

0.1254on79degreesoffreedomMultipleR-squared:

0.9935,AdjustedR-squared:

0.9931F-statistic:

2407on5and79DF,p-value:

2.2e-16上面语句第一行用来生成所需要的Q函数。

第二行用来生成线性回归模型中的时间输入，然后存储于叫做trend的数列中。

第三行语句是