R语言时间序列中文教程Word格式.docx
《R语言时间序列中文教程Word格式.docx》由会员分享，可在线阅读，更多相关《R语言时间序列中文教程Word格式.docx（43页珍藏版）》请在冰点文库上搜索。
1.运用R语言研究JJ数据
学习R言语时间序列分析程序操作，需要从最基础、最简单的学起，例如在命令窗口中，输入并执行2+2等于4的R语言命令。
2+2
[1]4
执行完2+2等于4的R语言命令后，我们可以开始时间序列，即学着把玩johnson&
Johnson数据。
下载jj.dat或执行下面语句。
这个数据已被人上传到因特网中。
R所需要做的只是将网址进行扫描就可以将数据读取进入R的编程环境中。
下面有3种不同读取数据的方法：
jj=scan（"
http:
//www.stat.pitt.edu/stoffer/tsa2/data/jj.dat"
）#readthedata读取数据
jj<
-scan（"
）#readthedataanotherway第二种方法读取数据
scan（"
）->
jj#andanother第三种方法读取数据
使用R语言的人，有的喜欢使用[<
-]，有的喜欢使用[->
]，大多数医疗系统的工作者喜欢用[=]，正因为如此才用了上面种不同读取数据的方法。
读取数据后，键入并执行jj，数据在窗口便会有如下显示：
jj
[1]0.710.630.850.44
[5]0.610.690.920.55
.....
[77]14.0412.9614.859.99
[81]16.2014.6716.0211.61
jj中有84个数据被称作对象。
下面命令可以显示所有对象。
objects（）
如果使用matlab，你会认为jj是一个84行1列的向量，但实际上不是这样。
jj有次序，有长度，但没维度，R称这些对象为向量，要小心区别。
在R里，矩阵有维度，但向量没维度。
这都是程序语言的一些概念。
jj[1]#thefirstelement列中第一个数据
[1]0.71
jj[84]#thelastelement列中最后一个数据
[1]11.61
jj[1:
4]#thefirst4elements列中第一至第四个数据
jj[-（1:
80）]#everythingEXCEPTthefirst80elements列中除第80个以外的所有数据
[1]16.2014.6716.0211.61
length（jj）#thenumberofelements有多少个数据
[1]84
dim（jj）#butnodimensions...但没维度
NULL
nrow（jj）#...norows没行
ncol（jj）#...andnocolumns没列
#如果你要把jj转变为一个向量，执行如下操作后，维度为84行1列。
jj=as.matrix（jj）
dim（jj）
[1]841
然后把jj转变为一个时间序列对象。
jj=ts（jj,start=1960,frequency=4）#ts（）命令
这个数据是从1960年开始，个个季度的收入，frequency=4指四个季度。
R语言的优势在于可用一条命令做很多事，即可以把前面的命令放在一起打包执行。
其操作如下：
jj=ts（scan（"
//www.stat.pitt.edu/stoffer/tsa2/data/jj.dat"
）,start=1960,frequency=4）
在上面命令里，scan可以被read.table替代。
用read.table读取数据可生成matrix对象，还可以给每列起名字。
下面学习一下read.table,dataframes,和时间序列对象。
输入命令后，窗口会有如下显示：
jj=ts（read.table（"
）,start=1960,frequency=4）
help（read.table）
help（ts）
help（data.frame）
需要注意的是，Scan和read.table不一样。
Scan生成的是有维度的向量，read.table生成的则是带有维度的数据架构。
读取jj数据的最后要领。
如果这个数据是从1960年第三个季度开始的，所需输入命令则为ts（x,start=c（1960,3）,frequency=4）；
如果是一个每月每月的数据，例如数据是从1984年6月开始的，需要输入的命令则为ts（x,start=c（1984,6）,frequency=12）。
输入命令后，转变后的时间序列对象为：
Qtr1Qtr2Qtr3Qtr4
19600.710.630.850.44
19610.610.690.920.55
.....
197914.0412.9614.859.99
198016.2014.6716.0211.61
注意到区别了吗？
时间信息，也就是4个不同的季度的数据被加载到里面了。
进行时间数据分析后，窗口会有如下显示：
time（jj）
19601960.001960.251960.501960.75
19611961.001961.251961.501961.75
......
19791979.001979.251979.501979.75
19801980.001980.251980.501980.75
接下来输入如下组合命令。
（jj=ts（scan（"
）,start=1960,frequency=4））
然后进行对数据绘图：
plot（jj,ylab="
EarningsperShare"
main="
J&
J"
）
输入以上命令后，可以看到如下结果：
再输入下面的命令，看看区别。
plot（jj,type="
o"
col="
blue"
lty="
dashed"
）
plot（diff（log（jj））,main="
loggedanddiffed"
下面利用操作plot.ts和ts.plot两个相关命令，显示区别。
x=-5:
5#sequenceofintegersfrom-5to5
y=5*cos（x）#guess
par（mfrow=c（3,2））#multifiguresetup:
3rows,2cols
#---plot:
plot（x,main="
plot（x）"
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:
ts.plot（x,main="
ts.plot（x）"
ts.plot（ts（x）,ts（y）,col=1:
2,main="
ts.plot（x,y）"
）#note-xandyaretsobjects
#---thehelpfiles[?
andhelp（）arethesame]:
?
plot.ts
help（ts.plot）
par#mightaswellskimthegraphicalparametershelpfilewhileyou'
rehere
从窗口中的显示可以看出，如果数据是时间序列对象，使用plot（）命令就足够了；
如果数据是平常序列，使用plot.ts（）也可以做时间绘图。
不过，把jj数据放在一张图上，数据会随着时间的变化上上下下跳动，能从整体上反应上升或者下降的趋势。
上文中用红色光滑的曲线代表上升的趋势，简单明了。
这需要将过滤和光滑的技巧使用在jj数据上。
在这里，我们用对称的移动平均值来达到过滤和光滑的目的。
下面使用公式：
fjj（t）=⅛jj（t-2）+¼
jj（t-1）+¼
jj（t）+¼
jj（t+1）+⅛jj（t+2）
除此之外，lowess的过滤平滑技巧（蓝色曲线）也要使用在jj数据中。
具体操作如下图：
k=c（.5,1,1,1,.5）#kisthevectorofweights用于对称移动平均的系数
（k=k/sum（k））
[1]0.1250.2500.2500.2500.125
fjj=filter（jj,sides=2,k）#?
filterforhelp[butyouknewthatalready]使用对称移动平均
plot（jj）
lines（fjj,col="
red"
）#addsalinetotheexistingplot称移动平均的绘图
lines（lowess（jj）,col="
）#lowess的绘图
操作后，窗口会显示下面结果：
看完jj数据，我们就需要开始具体分析。
第一步，我们把所有jj数据都取log值。
第二步，我们把log值做差，即使用log值数列中第二值减去第一值，第三值减去第二值，第四值减去第三值等等。
如果做差处理前数列里有n个数值，处理后的结果中将有n-1个数值。
dljj=diff（log（jj））#differencetheloggeddata做log和差的处理
plot（dljj）#plotitifyouhaven'
talready对结果制图
shapiro.test（dljj）#testfornormality测试结果的正态分布的性质
Shapiro-Wilknormalitytest
data:
dljj
W=0.9725,p-value=0.07211
处理完毕以上两步，我们接下来就要将柱形分布图和QQ图放在一起。
这两个图的本质仍旧在于测试数据正态分布的性质。
数据正态分布的性质是整个统计学构架坚实的基础，如果这个性质的存在比较可信、通过了所有测试，统计分析中得出的结论就比较可信、就通得过所有测试。
当然如果这个性质在数据中不存在，我们需要用其它的技巧来处理。
详细的，参看R语言样品比较应用的实例。
以上操作，在窗口中有如下显示：
par（mfrow=c（2,1））#setupthegraphics设置为两图的输出
hist（dljj,prob=TRUE,12）#histogram柱形分布图
lines（density（dljj））#smoothit-?
densityfordetails柱形分布图的曲线
qqnorm（dljj）#normalQ-QplotQQ图
qqline（dljj）#addaline在QQ图上加直线
经过测试数据后，窗口会有如下显示：
在实践操作中，时间序列数据存在着前后关系。
例如，今天股票的价格很有可能决定明天股票的价格。
明天的温度取决于今天的气温。
做天气预报的具体操作方法，是使用已经存在的天气历史记录，比如说今天的气温，昨天的气温，前天的气温等等，来预测明天的气温。
当然，在进行预测之前，我们一定要看清时间序列数据中的前后关系结构，清楚哪一个特定的历史数据可以精确预测未来的数据。
在这里，我们使用被log和求差后的dljj数据，来介绍分析时间序列数据前后关系结构的具体技巧。
在预测的实际应用中，我们总希望用历史数据来预测即将要产生的数据。
事实上，已产生的数据相对于即将产生的数据，中间存在着一定的延迟，也就是lag.比方说在某天凌晨12点开始记录室内温度，每小时记一次，一共连续记录了10个小时。
凌晨12点的数据和凌晨3点的数据之间就存在着延迟。
12点的数据比3点的早了3个小时，可记作lag3.3点的数据比12点的晚了3个小时，可记作lead3.
我们下面来介绍关联性。
例如，冷饮的销量与天气温度存在关联性。
温度越高冷饮销量就越高，温度越低冷饮销量也越低。
这种关联性称为正面关联性。
又如，人的体重和跑步速度也存在关联性。
不过，人的体重越重，跑步速度就会越慢，体重越低，相对来讲，速度就会越快。
这种关联性称为负面关联性。
下面我们回到预测应用上。
如果现在收集的数据与将来的数据之间存在着正面或者是负面的关联性，我们就可以用现在收集的数据来预测未来的数据。
因此找到现在收集到的数据与未来数据之间的关联性是最关键的。
找到这种关联性的具体技巧被称作延迟图表，也就是lag.plot.
我们可以在电脑里输入如下命令：
lag.plot（dljj,9,do.lines=FALSE）#whythedo.lines=FALSE?
...tryleavingitout
上面语句显示了lag.plot用法，输入的数据是被log和作差后的jj.dat数据。
其中9表示我们要考虑从1到9这9个不同的延迟。
值得注意的是，在下面图表中显示出延迟4和8显示出了正面关系。
其他几个延迟存在着负面关系。
我们可以利用这4和8的正面关系来进行预测，即用现有数据计算接下来的第4个或者是第8个数据。
下面我们来看ACF和PACF图表。
确定我们已经观察到的正面和负面关系。
par（mfrow=c（2,1））#Thepowerofaccurateobservationiscommonlycalledcynicism
#bythosewhohavenotgotit.-GeorgeBernardShaw
acf（dljj,20）#ACFtolag20-nographshown...keepreading
pacf（dljj,20）#PACFtolag20-nographshown...keepreading
#!
!
NOTE!
acf2onthelinebelowisNOTavailableinR...detailsfollowthegraphbelow
acf2（dljj）#thisiswhatyou'
llseebelow
在上图中，ACF和PACF横坐标的标记是1,2,3,4,5.但因为数据是季度性的，每年有4个季度所以1,2,3,4,5的标记代表的4,8,12,16,20的延迟。
当然，如果我们不喜欢上面横坐标的标记，也可以将dljj更改为ts（dljj,freq=1）;
也就是说acf（ts（dljj,freq=1）,20）。
因为在上面ACF图表中横坐标1代表的是延迟4，横坐标2代表的是延迟6，其相应的竖线代表的就是延迟4和8的正面关系。
接下来，下面我们介绍结构拆析。
在前面R代码中，我们曾将所有jj数据进行了lag变型。
在变型后的数据中，存在着上升趋势，季节的影响和每一时间点产生的随机的误差。
根据这一数据图，我们能够把趋势、季节和误差从变型后的jj数据中拆析出来，分别研究，或者分别进行绘图，以便于单独检查。
下面等式代表将要使用的数学模型：
Log（jj）=趋势+季节+误差
log（jj）=trend+season+error
结构拆析的R命令是stl（）,下面语句中stl命令中输入的是lag变型后的jj数据。
其中的“per”输入指的是使用季节循环来进行拆析。
stl语句在这里生成了一个叫dog的R物件，然后Plot语句将dog物件进行绘图。
具体操作如下图”
plot（dog<
-stl（log（jj）,"
per"
））
窗口会出现下面所示：
上图中有四行R的绘图，其中第一行代表原来的log（jj）的数据。
此数据可以看到总体的上升趋势还存在着一定季节循环性的变化。
第二行绘图代表拆析后季节循环的作用。
第三行绘图代表拆析后将季节循环作用清除剩余的上升趋势，此数据清楚地看到那种循环性变化已经不存在，剩余的只是趋势。
第四行绘图代表将季节循环作用和总体的趋势从数据中清除后所剩余的随机产生的误差。
如果我们需要对数据的误差进行一些常规检测，例如进行正态分布检测，绘制QQ图，还有绘制柱形图。
我们所需要的具体误差数据被存在叫做dog$time.series[,3]的数列里。
即叫dog的物件中有个叫time.series的数据矩阵，误差就被存储在这个数据矩阵的第三列里。
$指调取dog物件中的time.series数据矩阵。
[,3]指数据矩阵中第三列。
如果要对这一数列的误差值进行ACF的分析，只需要执行命令acf（dog$time.series[,3]）。
再接下来，我们对log变型后的jj数据进行线性回归模型分析。
与上面结构拆析不同的是，我们在这里使用四个季度来量化季节循环对数据的影响。
一年中有四个季度，也是我们所使用数据所代表的。
这个jj数据是某一家公司的季度收入数据，从上面绘图中我们就可以看到，每一年第三季度就会出现一个收入高峰，随之而来第四季度收入就会跌入低谷。
然后在一季度和二季度收入又会逐渐上升。
这也就是说，每一季度对这家公司收入的影响都是不一样的。
具体考虑到这种季度之间的不同，我们可使用如下数学模型：
log（jj）=β*time+α1*Q1+α2*Q2+α3*Q3+α4*Q4+ε
这个数学模型的意思是：
log（jj）=趋势*时间+一季度的影响*一季度+二季度的影响*二季度+三季度的影响*三季度+四季度的影响*四季度+误差
上面的模型β代表的就是总体上升趋势，α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（）用来输出模型参数数值。
具体操作以及显示如下：
Q=factor（rep（1:
4,21））#make（Q）uarterfactors[that'
srepeat1,2,3,4,21times]
trend=time（jj）-1970#notnecessaryto"
center"
time,buttheresultslooknicer
reg=lm（log（jj）~0+trend+Q,na.action=NULL）#runtheregressionwithoutanintercept
#--thena.actionstatementistoretaintimeseriesattributes
summary（reg）
Call:
lm（formula=log（jj）~0+trend+Q,na.action=NULL）
Residuals:
Min1QMedian3QMax
-0.29318-0.09062-0.011800.084600.27644
Coefficient
升级会员 