R语言时间序列中文教程.doc

39

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语言，还可以阅读优秀时间序列读物Ecomometrics in R，也可以上QuickR 网站。

目录

R语言时间序列分析 1
前言 1
目录 2
1.运用R语言研究JJ数据 3
2.运用R语言研究空气污染 19
3.运用R语言研究自动回归移动平均集成模型 22
4.运用R语言研究全球变暖理论 31
5.运用R语言研究非独立误差与线性回归 32
6.运用R语言研究估计波动的频率 36
7.运用R语言研究厄尔尼诺频率 39
8.运用R语言研究太阳黑子周期频率 40

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") # read the data读取数据
jj <- scan("http://www.stat.pitt.edu/stoffer/tsa2/data/jj.dat") # read the data another way第二种方法读取数据
scan("http://www.stat.pitt.edu/stoffer/tsa2/data/jj.dat") -> jj # and another第三种方法读取数据
使用R语言的人，有的喜欢使用[<-]，有的喜欢使用[->]，大多数医疗系统的工作者喜欢用[=]，正因为如此才用了上面种不同读取数据的方法。
读取数据后，键入并执行jj，数据在窗口便会有如下显示：
jj
[1] 0.71 0.63 0.85 0.44
[5] 0.61 0.69 0.92 0.55
. . . . .
. . . . .
[77] 14.04 12.96 14.85 9.99
[81] 16.20 14.67 16.02 11.61
jj中有84个数据被称作对象。下面命令可以显示所有对象。

objects()
如果使用matlab，你会认为jj是一个84行1列的向量，但实际上不是这样。jj有次序，有长度，但没维度，R称这些对象为向量，要小心区别。
在R里，矩阵有维度，但向量没维度。这都是程序语言的一些概念。
jj[1] # the first element列中第一个数据
[1] 0.71
jj[84] # the last element列中最后一个数据
[1] 11.61
jj[1:4] # the first 4 elements列中第一至第四个数据
[1] 0.71 0.63 0.85 0.44
jj[-(1:80)] # everything EXCEPT the first 80 elements列中除第80个以外的所有数据
[1] 16.20 14.67 16.02 11.61
length(jj) # the number of elements 有多少个数据
[1] 84
dim(jj) # but no dimensions ...但没维度
NULL
nrow(jj) # ... no rows 没行
NULL
ncol(jj) # ... and no columns没列
NULL
#如果你要把jj转变为一个向量，执行如下操作后，维度为84行1列。
jj = as.matrix(jj)
dim(jj)
[1] 84 1
然后把jj转变为一个时间序列对象。
jj = ts(jj, start=1960, frequency=4)#ts()命令
这个数据是从1960年开始，个个季度的收入，frequency=4指四个季度。
R语言的优势在于可用一条命令做很多事，即可以把前面的命令放在一起打包执行。其操作如下：
jj = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/jj.dat "), start=1960, frequency=4)
在上面命令里，scan可以被read.table替代。用read.table读取数据可生成matrix对象，还可以给每列起名字。
下面学习一下read.table, data frames, 和时间序列对象。输入命令后，窗口会有如下显示：
jj = ts(read.table("http://www.stat.pitt.edu/stoffer/tsa2/data/jj.dat"), 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)。
输入命令后，转变后的时间序列对象为：
jj
Qtr1 Qtr2 Qtr3 Qtr4
1960 0.71 0.63 0.85 0.44
1961 0.61 0.69 092 0.55
. . . . .
. . . . .
1979 14.04 12.96 14.85 9.99
1980 16.20 14.67 16.02 11.61
注意到区别了吗？时间信息，也就是4个不同的季度的数据被加载到里面了。
进行时间数据分析后，窗口会有如下显示：
time(jj)
Qtr1 Qtr2 Qtr3 Qtr4
1960 1960.00 1960.25 1960.50 1960.75
1961 1961.00 1961.25 1961.50 1961.75
. . . . . .
. . . . . .
1979 1979.00 1979.25 1979.50 1979.75
1980 1980.00 1980.25 1980.50 1980.75
接下来输入如下组合命令。
(jj = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/jj.dat"), start=1960, frequency=4))

然后进行对数据绘图：
plot(jj, ylab="Earnings per Share", main="J & J")
输入以上命令后，可以看到如下结果：

再输入下面的命令，看看区别。
plot(jj, type="o", col="blue", lty="dashed")
plot(diff(log(jj)), main="logged and diffed")
下面利用操作plot.ts和ts.plot两个相关命令，显示区别。
x = -5:5 # sequence of integers from -5 to 5
y = 5*cos(x) # guess
par(mfrow=c(3,2)) # multifigure setup: 3 rows, 2 cols
#--- 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.plotx)")
ts.plot(ts(x), ts(y), col=1:2, main="ts.plot(x,y)") # note- x and y are ts objects
#--- the help files [? and help() are the same]:
?plot.ts
help(ts.plot)
?par # might as well skim the graphical parameters help file while you're here

从窗口中的显示可以看出，如果数据是时间序列对象，使用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) # k is the vector of weights用于对称移动平均的系数
(k = k/sum(k))
[1] 0.125 0.250 0.250 0.250 0.125
fjj = filter(jj, sides=2, k) # ?filter for help [but you knew that already]使用对称移动平均
plot(jj)
lines(fjj, col="red") # adds a line to the existing plot称移动平均的绘图
lines(lowess(jj) col="blue", lty="dashed")#lowess 的绘图
操作后，窗口会显示下面结果：

看完jj数据，我们就需要开始具体分析。
第一步，我们把所有jj数据都取log值。第二步，我们把log值做差，即使用log值数列中第二值减去第一值，第三值减去第二值，第四值减去第三值等等。如果做差处理前数列里有n个数值，处理后的结果中将有n-1个数值。
dljj = diff(log(jj)) # difference the logged data做log和差的处理
plot(dljj) # plot it if you haven't already对结果制图
shapiro.test(dljj) # test for normality 测试结果的正态分布的性质
Shapiro-Wilk normality test
data: dljj
W = 0.9725, p-value = 0.07211
处理完毕以上两步，我们接下来就要将柱形分布图和QQ图放在一起。这两个图的本质仍旧在于测试数据正态分布的性质。数据正态分布的性质是整个统计学构架坚实的基础，如果这个性质的存在比较可信、通过了所有测试，统计分析中得出的结论就比较可信、就通得过所有测试。当然如果这个性质在数据中不存在，我们需要用其它的技巧来处理。详细的，参看R语言样品比较应用的实例。
以上操作，在窗口中有如下显示：
par(mfrow=c(2,1)) # set up the graphics 设置为两图的输出
hist(dljj, prob=TRUE, 12) # histogram柱形分布图
lines(density(dljj)) # smooth it - ?density for details柱形分布图的曲线
qqnorm(dljj) # normal Q-Q plot QQ图
qqline(dljj) # add a line 在 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) # why the do.lines=FALSE? ... try leaving it out
上面语句显示了lag.plot用法，输入的数据是被log和作差后的jj.dat数据。其中9表示我们要考虑从1到9这9个不同的延迟。值得注意的是，在下面图表中显示出延迟4和8显示出了正面关系。其他几个延迟存在着负面关系。
我们可以利用这4和8的正面关系来进行预测，即用现有数据计算接下来的第4个或者是第8个数据。

下面我们来看ACF和PACF图表。 确定我们已经观察到的正面和负面关系。
par(mfrow=c(2,1)) # The power of accurate observation is commonly called cynicism
# by those who have not got it. - George Bernard Shaw
acf(dljj, 20) # ACF to lag 20 - no graph shown... keep reading
pacf(dljj, 20) # PACF to lag 20 - no graph shown... keep reading
# !!NOTE!! acf2 on the line below is NOT available in R... details follow the graph below
acf2(dljj) # this is what you'll see below

在上图中，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代表的是四个季度的影响。有一个非常有趣的问题，上面模型是把所有四个季度的趋势都加在了一起，其结果却是某单一季度的收入。四个季度的和如何能够与一个季度相等问题就出在Q1 Q2 Q3 Q4 上。因为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)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)
Call:
lm(formula = log(jj) ~ 0 + trend + Q, na.action = NULL)

Residuals:
Min 1Q Median 3Q Max
-0.29318 -0.09062 -0.01180 0.08460 0.27644

Coefficients:
Estimate Std. Error t value Pr(>|t|)
trend 0.167172 0.002259 74.00 <2e-16 ***
Q1 1.052793 0.027359 38.48 <2e-16 ***
Q2 1.080916 0.027365 39.50 <2e-16 ***
Q3 1.151024 0.027383 42.03 <2e-16 ***
Q4 0.882266 0.027412 32.19 <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
上面语句第一行用来生成所需要的Q函数。第二行用来生成线性回归模型中的时间输入，然后存储于叫做trend的数列中。第三行语句是用lm()命令来建立线性模型。在R语言中，~用来分隔进行预测的变量和被预测的变量。
在上面语句第三行中，log（jj）是被预测的变量，它被放置在~的左边。而且0，trend, 和Q函数都是用来进行预测的变量。它们被放在~右边。其目的就是要用0TREND和Q function来预测Log（jj）。
现在，我们复习一下什么叫做线性。线性是在坐标平面上的一条直线。坐标平面上有两个数轴，横轴叫做x轴，竖轴叫做y轴。在平面上画一条直，可使用截距和斜率定义它。如果这条直线通过原点，其截距为0. 在~右边的0是我们模型中的截距，也就是说我们模型中的直线通过原点。如果我们决定使用有截距的线性模型，即模型中的直线不通过原点，就要把第三条R语句中的0改为1.这样，R语言就可以将截距计算出来。
执行summary（）命令后，我们就可得到具体模型参数的数值。模型中的趋势，即β估计值为0.167172，一季度的影响记为a1，估计值为1.052793，第二季度的影响记为a2，估计值为1.080916，第三季度的影响记为a3，估计值为 1.151024，第四季度的影响记为a4，估计值为 0.882266。
随后，我们将所有估计值插入上面模型中，可得到下面等式：
Log(jj)=0.167172*时间+1.052793*Q1+ 1.080916*Q2+ 1.151024 *Q3+0.882266*Q4+误差
然后，我们就可以进行预测比如说，在第200个时间所处的季度为4，其预测值就可被计算为：
Log(jj)= 0.167172*200+1.052793*0+ 1.080916*0+ 1.151024 *0 +0.882266*1=Now check out what happened. Look at a plot of the observations and their fitted values:
在上面，我们已经知道如何使用计算出来的模型对某一时间点进行预测。但是，我们的预测是否符合在实际中观察到的时间序列数据，我们需要将实际当中观察到的数据与模型计算出来的数据放在一起，绘成图表进行比较。
接下来要进行的工作是，将预测的数据和实观察的数据进行比较。下面第一条语句是对实际观察到的log(jj)数据进行绘图。第二行语句是对计算出的数据进行绘图，其中fitted(reg)是使用我们已建立好的模型来计算预期值。
plot(log(jj), type="o") # the data in black with little dots
lines(fitted(reg), col=2) # the fitted values in bloody red - or use lines(reg$fitted, col=2)
输入上面两个命令后，窗口就会显示如下图表。在图表中，黑线代表实际中观察到的log(jj)数据。红线代表由模型计算出来的log(jj)数据。它们看上去非常相似，但区别在于模型中计算出来的数值没误差。

... and a plot of the residuals and the ACF of the residuals:
在对log(jj)建立模型时，具体的误差也被估计了出来，并存储于一个叫做residuals的数列中，可使用命令resid(reg)进行调用。
下面，plot(resid(reg))语句是用误差的绘图来确定误差的变化是否比较小范围之内。acf(resid(reg),20)语句是用来检验误差的自动关联性。因为我们知道误差之间不应当存在任何关联性。如果误差之间没关联性存在，只有在lag0的ACF上有数值为1的关联性。具体操作和显示如下：
par(mfrow=c(2,1))
plot(resid(reg)) # residuals - reg$resid is same as resid(reg)
acf(resid(reg),20) # acf of the resids
下面给出了相关误差的两个绘图。第一个图可以看到误差的变化范围在允许范围之内，并且没大变化。第二张图只有在lag0上有数值为1的关联性，其他lag上的关联性都非常小，即说明误差和误差之间是没关系。这符合我们模型所需要假设的。

2.运用R语言研究空气污染

接下来，我们讲解一个延迟线性回归模型，同样要使用到R语言中lm()命令来生成模型。我们以空气污染研究为实例进行讲解。
在实例中，我们要使用空气中的固体污染物含量来预测心血管疾病死亡人数，并且延迟四周，即一个月。心血管疾病死亡人数的数据记为cmort.dat, 空气中固体污染物的数据记为part.dat。先将相关数据传到网上，然后用R语言扫描其连接，进而读取数据。
下面是R语言读取两个数据所使用的命令语言。
mort = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/cmort.dat"),start=1970, frequency=52) # make these time series objects
Read 508 items
part = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/part.dat"),start=1970, frequency=52)
Read 508 items
读取数据后，运用ts.intersect()命令把mort 数据、part数据以及part经过延迟四周的数据捆绑在了一起，使三列数据形成行与行之间的正确对应关系，然后进行接线性回归分析。
在这个线性模型中，我们要使用当前空气污染物的数据和四周以前空气污染物数据对当前心血管病死亡的人数进行预测。因为我们使用四个星期以前的数据，因而所建立起来的模型被称为延迟模型，由R语言中的lag（）命令来实现，延迟四个星期则用-4表示。
ded = ts.intersect(mort,part,part4=lag(part,-4), dframe=TRUE) # tie them together in a data frame
接下来的语句表示的是如何生成线性模型，使用的也是lm（）命令。在语句中，mort数列表示被预测的变量，放在~左边，part变量和延迟了四周的part变量属于用来进行预测的变量，要放在~右边。上面已被捆绑好的数据被称作ded，要跟在data = 后面。
fit = lm(mort~part+part4, data=ded, na.action=NULL) # now the regression will work
接下来使用summary命令输出所有估计完成的参数。
summary(fit)
Call:
lm(formula = mort ~ part + part4, data = ded, na.action = NULL)

Residuals:
Min 1Q Median 3Q Max
-22.7429 -5.3677 -0.4136 5.2694 37.8539

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 69.01020 1.37498 50.190 < 2e-16 ***
part 0.15140 0.02898 5.225 2.56e-07 ***
part4 0.26297 0.02899 9.071 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.323 on 501 degrees of freedom
Multiple R-Squared: 0.3091,