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语言，还可以阅读优秀时间序列读物Ecomometrics in R，也可以上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种不同读取数据的方法:jj scan " # read the data读取数据jj - scan " # read the data another way 第二种方法读取数据scan " - 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.71jj[84] # the last element列中最后一个数据 [1] 11.61jj[1:4] # the first 4 elements列中第一至第四个数据 [1] 0.71 0.63 0.85 0.44jj[- 1:80 ] # everything EXCEPT the first 80 elements列中除第80个以外的所 有数据 [1] 16.20 14.67 16.02 11.61length jj # the number of elements 有多少个数据 [1] 84dim jj # but no dimensions ...但没维度 NULLnrow jj # ... no rows 没行 NULLncol 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 " " , start 1960, frequency 4 在上面命令里，scan可以被read.table替代。用read.table 读取数据可生成matrix对象，还可以给每列起名字。 下面学习一下read.table, data frames, 和时间序列对象。输入命令后，窗口会有如下显示: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 。 输入命令后，转变后的时间序列对象为:jj Qtr1 Qtr2 Qtr3 Qtr4 1960 0.71 0.63 0.85 0.44 1961 0.61 0.69 0.92 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 " , 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" 下面利用操作两个相关命令，显示区别。x -5:5 # sequence of integers from -5 to 5y 5*cos x # guesspar 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.plot x " 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.tshelp 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.125fjj 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处理完毕以上两步，我们接下来就要将柱形分布图 和图放在一起。这两个图的本质仍旧在于测试数据正态分布的性质。数据正态分 布的性质是整个统计学构架坚实的基础，如果这个性质的存在比较可信、通过了 所有测试，统计分析中得出的结论就比较可信、就通得过所有测试。当然如果这 个性质在数据中不存在，我们需要用其它的技巧来处理。详细的，参看R语言样 品比较应用的实例。以上操作，在窗口中有如下显示:par mfrow c 2,1 # set up the graphics 设置为两图的输出hist dljj, prob TRUE, 12 # histogram柱形分布图 lines density dljj # smooth it - ?density for details柱形分布图的曲线 norm dljj # normal Q-Q plot 图 line dljj # add a line在 图 上 加 直线 经过测试数据后， 窗口会有如下显示:在实践操作中，时间序列数据存在着前后关系。例如，今天 股票的价格很有可能决定明天股票的价格。明天的温度取决于今天的气温。做天 气预报的具体操作方法，是使用已经存在的天气历史，比如说今天的气温， 昨天的气温，前天的气温等等，来预测明天的气温。当然，在进行预测之前，我 们一定要看清时间序列数据中的前后关系结构，清楚哪一个特定的历史数据可以 精确预测未来的数据。在这里，我们使用被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 Shawacf dljj, 20 # ACF to lag 20 - no graph shown... keep readingpacf 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 belowacf2 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)的数据。此数据可以看到总体的上升趋势还存在着一定季节循环性的变化。第二行绘图代表拆析后季节循环的作用。第三行绘图代表拆析后将季节循环作用清除剩余的上升趋势，此数据清楚地看到那种循环性变化已经不存在，剩余的只是趋势。第四行绘图代表将季节循环作用和总体的趋势从数据中清除后所剩余的随机产生的误差。 如果我们需要对数据的误差进行一些常规检测，例如进行正态分布检测，绘制图，还有绘制柱形图。我们所需要的具体误差数据被存 在叫做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 nicerreg lm log jj ~0+trend+Q, na.action NULL # run the regression without an intercept#-- the na.action statement is to retain time series attributessummary reg Call:lm formula log jj ~ 0 + trend + Q, na.action NULL Residuals: Min 1Q Median 3Q -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 freedomMultiple 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 " ,start 1970, frequency 52 # make these time series objects Read 508 items part ts scan " ,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 -22.7429 -5.3677 -0.4136 5.2694 37.8539 Coefficients: Estimate Std. Error t value Pr |t| Intercept 8 50.190 2e-16 *** part 0.15140 0.02898 5.225 2.56e-07 *** part4 0.26297 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, Adjusted R-squared: 0.3063 F-statistic: 112.1 on 2 and 501 DF, p-value: 2.2e-16 如上图所示，我们便获得了所研究需 要的数据。3.运用R语言研究自动回归移动平均集成模型 现在来介绍自动回归 移动平均集成模型，即ARIMA模型。AR指自动回归，MA指移动平均，I指集成。 我们先来生成两个随机的自动回归数列，记为x1和x2。它们都是ar1数列，因 为生成过程中使用order c 1,0,0 的参数输入。X1的AR参数是+0.9，x2的AR 参数是-0.9。 两个语句生成完毕后，我们以将这两个数列绘图出来。两个时间 序列非常不同，x1的参数为+0.9，它的图像类似于股票价格的变动图形;x2的 参数为-0.9， 它的图像看上去像是一列声波。窗口具体显示如下: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 2,1 plot x1 plot x2 在R语言中， ACF指自动关系性，ACF即Auto-Correlation Function的简称. 比方说，股票 价格今天的价格跟昨天的价格有关系，明天的价格会跟今天的或者昨天的价格有 关系。它们之间的关系性便用ACF来衡量。PACF被称作不完全自动关系性。自 动关系性ACF中存在着线性关系性和非线性关系性。不完全自动关系性就是把线 性关系性从自动关系性中消除。如果在线性关系性被去除以后，两个时间点之间 的关系性也就是不完全关系性。当PACF近似于0，这表明两个时间点之间的关 系性是完全由线性关系性所造成的。如果不完全关系性在两个时间点之间不近似 于0，这表明线性模型是不能够表达这两个时间点之间的关系。我们用下面的语句来检验x1和x2数列中的ACF和PACF。par mfcol c 2,2 acf x1, 20 acf x2, 20 pacf x1, 20 pacf x2, 20 对数列x1来讲，ACF从lag1到lag12都高于蓝色虚线，也就是说两个时间点的距离在1到12之间它们的自动关系都是正面的。所有线性关系在x1数列中被清除了结果是Partial ACF的x1图表。在PACF x1图表上可以看到只有lag0值为1，其他的lag上的关系值都低于蓝色虚线，近似于0，也就是说在x1数列中存在的自动关系基本上是线性关系。使用线性模型可以符合x1数列的要求。对数列x2来讲，两个时间点之间的ACF关系有负面的也有正面的，但经过去除线性关系以后所有不完全自动关系都接近于0，也就是说x2数列中的自动关系也基本上是线性的。使用线性模型符合x2数列要求。下面，我们来介绍移动平均MA模型。我们先用ARIMA.SIM 命令生成一个移动平均的数据列。不同的是，order c 0,0,1 ，ma1代表生成的数列，ma是一个参数，设置为+0.8，数列中有100个数据。其余3行命令是用来设置输出窗口绘制图表和计算ACF，PACF数值的。 经过一系列操作后，代码及其结果如下所示:# an MA1 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 acf x,20 pacf x,20 需要注意的是，因数据都是随机生成的，我们在重新运行这些命令时，生成的数据会与以前所生成的数据完全不同，其结果也不一样。除非在以上代码的前面使用随机锁定的命令set.seed 1 ，每次生成的就是相同的数据了。当然也可以用set.seed 100 ，和其它数字输入都是可以的。在上图中，MA1代表了数列输出的结果。最上面显示的是数列的数目，横向坐标从1到100代表有100个数据。中间的图显示的是x数列的自动关系性，在lag0上有很强的自动关系性，因为Lag0表示没 时间延迟，当前数据与它有完全的关系性。所以lag0上的关系性什么都不代表。 而Lag1上存在着正面的关系性，即当前数据与前一个数据的关系性是正面的。 其他延迟上的自动关系性都近似于0，可以不被考虑。 order c 2,0,0 代表AR2 模型。下面，我们来看生成AR2数列及分析结果。在这里，需要有两个AR参数 1和-0.9。输入相关命令后，小窗会有如下显示:# an AR2 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 接 下来，我们介绍一下如何生成ARIMA(1,1,1)模型的随机数列。其中的设置order c 1,1,1 说明该数列是ARIMA模型，且AR的参数为0.9，MA的参数为0.5，数 列长度为200。输入相关命令后，小窗会有如下显示:# an ARIMA 1,1,1 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 我们再介绍一下自动回归移动平均集成模型的估计分析，即 ARIMA模型。也就是说，我们手头上有一个时间序列数据，如果知道这个数据是 被ARIMA模型生成的，我们希望能估计具体的ARIMA模型的参数，可以对整个数 据生成的过程进行重演和预测。下面，我们随机生成一些ARIMA模型的数据，并 且假装我们不知道模型的参数，然后做估计练习。 生成ARIMA数据的命令为 arima.sim . 自动回归的参数为0.9，移动平均的参数为-0.5，估计中要假设 这两个参数的设置为未知，并进行估计。计算估计值的命令为arima 。其中， x指随机生成的数据。在下面的输出中，我们可以看到AR参数的估计值为0.8465， ma的参数估计值为-0.5021，这两个估计值与前面生成数据时设置的参数非常相 似。也就是说，估计得比较精确。具体操作和显示如下所示: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 resultsCall: arima x x, order c 1, 0, 1 Coefficients: ar1ma1 intercept -- NOT the intercept - see Issue 1 0.8465 -0.5021 0.5006 s.e. 0.0837 0.1356 0.3150 sigma^2 estimated as 1.027: log likelihood -143.44, aic 294.89 下面来介绍衡量参数估计的可信度的方法。我们对估计的模型进 行检测和评估，需要执行的命令为tsdiag 。其中，x.fit指代估计出来的模型。 这个语句输出的结果如以下几个图表所示。第一个图表代表估计模型误差的绘 图。英文叫做Standardized Residuals, 上面有很多竖线在横向坐标的上下分 布。如果这个估计的模型比较可信，竖线的长度是比较相似的。如果竖线的长度 互相有很大出入或者根本就不同，估计模型的可信度就非常差。下面误差绘图中 竖线的长度比较相似，都处在稳定范围之内，即估计的模型没产生不符合要求的 误差分布。再介绍输出的第二张绘图，标题是ACF of Residuals。ACF指数据点 相互之间的关系，当然在生成这个数据时，数据点之间互相独立，并不存在任何 关系。所以在这张图上，只有位于0刻度上的竖线最高，其ACF值为1。 这个0 代表数据点与自己相比较， 即数据点永远和它自己有关系，这种关系数值为1。 其他横向数轴上的刻度代表一个数据点于其他数据点之间的关系，这些刻度上竖 线的长度几乎等于0，即这个数据点与其他数据点没明显关系。这张ACF图代表 估计的模型没造成误差之间的任何关系。这是符合数据生成时每个数据都是独立 的这个前提的。由此可见，这ACF图符合检测要求。下面来介绍第三张图，也就 是Ljung-Box 指标。这个指标可对每一个时间序列的延迟进行显著性的评估。 这张图的横坐标代表时间序列的延迟，纵坐标代表P-value，即显著性。如果 P-value十分小，就说明在其相对应的延迟点上是显著的。我们就需要抛弃所假 设的模型，并且结论在所假设的模型不可信。需要注意的是，他们使用假设的模 型对一个时间序列进行估计，如果P-value是显著的话，我们使用的模型就不可 信，需要尝试其他新模型。具体判定技巧是，P-value点的高度越高，我们的模 型越可信。 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在对模型的参数进行了估计后，我们接下来要 对这些估计的参数进行预测。我们需要做的是向前预测10个数据点。预测所使 用的语句是predict . 在 中输入x.fit 和n.ahead 10 ，其中x.fit 代表 前面估计出来的模型与参数，n.ahead 10 代表向前预测10个数据点。这条语 句给我们生成的预测程序项目，被命名为x.fore. 在下面两条语句中，计算了 两列的新数据，分别被称为U和L. U指最高界限，L指最低界限。 最高限和最 低限告诉我们未来数据有很高可能性事发生在两个界限之间，超出或低于这两个 界限的几率并不高。 下一步使用Min 和 条语句，把数据中的最大值和最小值 挑选出来，为接下来的统计绘图设置图标界限。再使用TS.PLOT 语句. 把原有 的数据和向前预测10个数据点的数据绘在一张图上。最后使用lines 语句， 把计算的最高界限和最低界限加到绘好的图上。这两条界限的颜色为蓝色，形式 也被设置为虚线的形式。x.fore predict x.fit, n.ahead 10 # plot the forecastsU x.fore$pred + 2*x.fore$seL x.fore$pred - 2*x.fore$seminx min x,L x x,U ts.plot x,x.fore$pred,col 1:2, ylim c minx,x lines U, col "blue", lty "dashed" lines L, col "blue", lty "dashed" 下面的绘图是 由上面的语句所完成的。图右侧除去原有的数据外，又加上了三条曲线。最上面 和最下面的曲线是虚线，是计算出的最高界限和最低界限。中间的曲线为实线， 是对未来的10个数据点进行的预测。在计算最高、最低限时，我们使用了一个 2的系数，也就是说，这两个界限之间的几率为95%。 在未来的10个数据点有 95%的几率是出现在上限和下限之间的。 而出现在两个界限之外的几率小于或等 于5%。 4.运用R语言研究全球变暖理论下面，来介绍R语言在全球温度变暖研 究中的应用。我们将全球温度变暖研究的数据叫globtemp2.dat. 数据中有3列， 第二列指全球每年温度从1880年开始到2004年。这个数据已经在网上已公布， 我们只需要用下面语句读取就行。第一条语句读取全球温度数据生成一个叫做u 的数列矩阵;第二条语句读取的第二列数据变成时间序列，命名为gtemp;第三 条语句是将gtemp数据绘成图来看温度变化的具体情况。操作以后，窗口里会有 如下显示:u read.table " # read the datagtemp ts u[,2], start 1880, freq 1 # yearly temp in col 2plot gtemp # graph the series not shown here 经过多年，科学家们论证这个时间序列符合ARIMA 1,1,1 模型， 并有0.6 oC每世纪的漂移。这也就是所谓全球变暖理论。下面我们要计算出具 体模型。我们有近200年的气温数据，且知道具体模型是ARIMA，下面我们需要 估计AR和MA的具体数值。语句如下，我们从结果中可以看到AR的数值为0.2545， MA的数值为-命令具体显示如下:arima gtemp, order c 1,1,1 Coefficients: ar1 ma1 0.2545 -0.7742 s.e. 0.1141 0.0651 下面，我们就来介绍 漂移参数的估计方法。具体操作如下所示:drift 1:length gtemp arima gtemp, order c 1,1,1 , xreg drift Coefficients: ar1 ma1 drift 0.2695 -0.8180 0.0061 s.e. 0.0030在上面的语句中，第一行生成了一个名为drift 的空数列，其长度与全球温度数列的长度相同。而在第二个语句中，我们只要把 drift数列作为xreg输入就可以了。从输出中，我们可以看到AR参数值为 0.2695，MA参数值为-0.8180，是drift值为0.0061，即漂移产生的值. 上面提 到过漂移值为+0.6，我们的估计值比提供的数值小100倍。其原因是因为我们使 用的全球温度数据是经过被100整除的，所以估计的漂移值也相应地小了100 倍。如果将gtemp数列输出后看一下，我们就可以明白所有温度都是被100整除 过的，就不会有什么疑惑了。具体情况，代码显示如下:u read.table " # read the datagtemp ts u[,2], start 1880, freq 1 # yearly temp in col 2plot gtemp drift 1:length gtemp arima gtemp, order c 1,1,1 , xreg drift gtemp5.运用R语言研究非独立误差与线性回归 在前面我们介绍的模型当 中，误差都是相互独立的，即误差之间不互相影响。但如果模型中存在着互相影 响的误差，我们运用R语言该怎么做呢,在下面例子当中，误差是一个自动回归 2的数列，记作AR2，意思是现在的误差依赖于前两个已经被观测到的误差,误差 之间是非独立的。具体有两个意思，我们要学习如何从数据中确认所有的误差符 合自动回归模型，而且参数为2;我们要学习如何把这个AR2模型的误差加入到 线性回归模型当中。我们使用的数据仍旧是空气中污染物和心血管疾病死亡率的 数据。在下面模型中，我们认为，时间可以决定心血管疾病死亡人数，随着时间 向前移动，心血管疾病死亡人数会呈现下降趋势。对于空气中的污染物来讲，空 气中的污染物越高，心血管疾病的死亡人数也会越高。所以在这个模型中，心血 管疾病死亡人数是被预测的变量。时间变量和空气中污染物的变量是用来进行预测的变量。 第一句命令是调用nlme，R语言程序包。其具体显示如下: library nlme # load the package第二句命令是将时间数列从mort数列当中提出，用来作为进行预测的变量称之为trend。命令如下: trend time mort # assumes mort and part are there from previous examples 第三句命令，先建立线性回归模型，再从模型中提取模型的误差来研究误差之间的相互关系性，以便确认误差之间是相互独立的，还是符合前面提到的AR2模型。因为模型建立以后，我们并不知道模型的误差是否是相互独立的。fit.lm lm mort~trend + part # ols第四行命令使用上面第三行命令中生成的被称作 fit.lm线性回归模型结果。 命令Resid 是用来提取线性回归模型中的误差数列的。Acf()命令是用来研究线性模型误差的自动关系性的，也就是确定误差之间是否相互关系。 acf resid fit.lm # check acf and pacf of the resids第五行命令同样是将模型中误差数列提出后使用PACF 命令进行不完全相关性的研究。 pacf resid fit.lm # or use acf2 resid fit.lm if you have acf2上面的pacf绘图是第五行命令运行的结果。在上图中，横向坐标代表的延迟，被叫做lag. 在lag1和lag2上出现了高于蓝色虚线的pacf,其他lag上超出蓝色虚线的距离并不明显可以忽略。这两个高于蓝色虚线的pacf说明这个误差数列符合AR2模型。第六行语句，使用gls 命令来对数据建立带有AR2误差的线性回归模型，生成结果记为fit.gls。 命令gls 在这里使用时间和空气污染物来预测心血管疾病死亡人数，然corARMA(p 2)是告诉R语言模型的误差符合AR2模型，如果模型的误差符合MA2模型就要改为corARMA q 2 . 其中method ”ML”,代表使用极大似然方法进行模型估计。大似然方法是一个 统计学的核心理念，极为重要，需要使用相关的教材着重学习。具体操作如下: # resids appear to be AR 2 ... now use gls from nlme:fit.gls gls mort~trend + part, correlation corARMA p 2 , method "ML" 第七行语句使 用summary命令来输出上面fit.gls模型结果的具体参数值。从下面窗口显示中， 我们可以看到这个线性模型的截距为3131.5452，模型中时间对于心血管疾病死 亡人数的影响，即时间的斜率值为-1.5444，也就是说随着时间每增长一周，心 血管疾病死亡率都要降低1.5444。 模型中空气污染物对心血管疾病死亡人数的 影响为0.1503，也就是说空气中污染物每上升一个单位，心血管疾病死亡人数 会增加0.1503个人。下面是电脑里显示的内容:summary fit.gls Parameter estimate s : Phi1 Phi2 0.3980566 0.4134305 Coefficients:Value Std.Error t-value p-value Intercept 3.653166 3e-04 trend -1.5444 0.4340 -3.558021 4e-04 part 0.1503 0.0210 7.162408 0e+00 # resid analysis- we assumed et φ1 et-1 + φ2 et-2 + wt where wt is weight. 事实上，ARIMA命令也可以完成上述自动关系误差线性模型。使用ARIMA命令， 其中输入的第一个参数是mort数列，表示Mort 是被预测的变量，第二个输入 的是Order c 2,0,0 代表AR2误差模型，xreg cbind trend,part 代表用时间 和空气污染物作为进行预测的变量。 下面我们用ARIMA命令从新估计这个模型， 所输出的模型参数可以看出截距为3132.7085，时间的影响是-1.5449，空气污 染物的影响是0.1503。这些数值与我们上面得到的结果之间略有差异是由计算 中四舍五入造成的。 fit2.gls arima mort, order c 2,0,0 , xreg cbind trend, part Coefficients: ar1ar2 intercept trend part 0.3980 0.4135 3132.7085 -1.5449 0.1503 s.e. 0.0405 0.0404 854.6662 0.4328 0.0211 sigma^2 estimated as 28.99: log likelihood -1576.56, aic 3165.13 Box.test resid fit2.gls , 12, type "Ljung" # and so on ... 6.运用R语言研究估计波动的频率我们来介绍如何运用R语言来估计波动的频 率。波动是一种非常常见的物理现象。人说话发出声波，手机通讯通过电波。我 们可以用R语言对一段被记录下来的波动进行频率的估计。比方说，每一年都会 有流感出现，得流感的人数在一年中的某个时间都是不一样的。对流感人数的记 录，就是对流感周期的记录。我们能够知道流感发生的频率，有利疾病的防治。 又如，股票市场其价格也是存在周期的，如果我们能计算出股票价格变化的周期， 就可以预测市场的变化。我们随机生成一个波动的时间序列数据来进行分析，阐 述R语言在估计拨动的频率。 下面的语句是用来生成一些自动回归2，即AR2 波动数据的。具体操作如下: x arima.sim list order c 2,0,0 , ar c 1,-.9 , n 2^8 # some data第二行R语句是对设定的自动回归模型的复杂数根进行计 算。这两个复数根的长度都小于1，值得注意，这是自动回归模型特有的性质。 具体操作如下所示: u polyroot c 1,-1,.9 # x is AR 2 w/complex roots 因为我们已经知道数据是从自动回归模型中生成的，根据模 型条件，对于其频率可进行精确的计算，并不需要估计。下面语句就是对这个波 动频率的计算，其模型已知为自动回顾模型。经过四舍五入，频率 u[1] / 2*pi # dominant frequency around .16: 下面这个语句是对绘图输出窗口的一些具 体设置。如果不进行这个设置，在绘图输出窗口中，每次只能输出一张图。我们 需要在窗口中一次输出3张图，并且罗列在一起。这3张图罗列在一起就形成了 3行1列的布局，也就是下面语句中要使用 3,1 的输入。 par mfcol c 3,1 在 估计频率之前，需要看一看随机生成数据的样子。需要把这个数据绘到图上。生 成绘图的语句如下: plot.ts x 接下来，我们对已生成的随机波动数据的频率 进行估计。注意，在这里不能够使用自动回归模型作为已知条件来进行频率的计 算，而是要对频率进行估计。这里使用的R语言是SPEC.PGRAM .在下面语句中， X代表生成的随机波动数据，SPANS C(3,3)是对计算的侧重进行的设置。 spec.pgram x, spans c 3,3 , log "no" # nonparametric spectral estimate; also see spectrum 如果使用者对上面的语句有疑问，可将其介绍文件查调出 来。具体使用语句如下:?spec.pgram # some help 在生成随机波动数据时，我们对于自动回归模型的输入设置为2, 即AR 2 。假 设我们已知模型为自动回归，但不知道其设置为2，我们也可以使用不完全的信 息进行频率估计。但是由于不知道其设置为2，所以频率也只能是估计出来的。 具体输入命令为:spec.ar x, log "no" # parametric spectral estimate下面是执行上面语句后的输出结果。我们来看窗口中输出的绘成图3 个图表。最上面图表显示的是随机波动的走势。中间的图表是在假设，用自动回 归模型的条件未知下进行的估计。在最顶端峰相对应的水平数轴的值就是我们所 估计出来的频率。具体来讲，水平数轴上标有刻度，.山峰最顶端所对应的刻度 在0.16左右，与前面所计算出的理论值相对应。最下面的绘图与中间的绘图结 论相似。峰值所对应的水平轴刻度也是在0.16左右。但这个绘图是在已知自动 回归模型，未知模型设置为2的条件下估计出来的。所以与中间的估计略有不同， 只是更精确。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 [1] 0.5555556+0.8958064i 0.5555556-0.8958064i Arg u[1] / 2*pi # dominant frequency around .16: [1] 0.1616497 par mfcol c 3,1 plot.ts x spec.pgram x, spans c 3,3 , log "no" # nonparametric spectral estimate; also see spectrum ?spec.pgram # some help spec.ar x, log "no" # parametric spectral estimate输出的图表:7.运用R语言研究厄尔尼诺频率下面，我们介绍用R语言研究厄尔 尼诺现象的规律。从1950年到1987年，每一个月都会在中太平洋地区收集空气 压力的变化，用来研究海水表面温度的变化。所收集的数据被称为南部波动指数， 英文叫做Southern Oscillation Index SOI ，因太平洋海水的温度每3到7年 就变暖一次，人称作厄尔尼诺现象。我们运用R语言命令直接读取数据，并依次 利用下面的语句研究厄尔尼诺现象变化的周期。 第一行R语句，是用来读取SOI 数据的。具体情况如下图: soi scan " 第二行R语句，是用来读取另外一 个数据的，被称作是rec。rec是中部太平洋鱼类数目每月的数据。rec的月份 与SOI数据的月份相对应。具体操作如下: rec scan " 第三行语句，是对 SOI数据的变化频率进行估计。操作如下:soi.per spec.pgram soi, taper 0, log "no" 第四条和第五条语句，是在1/12和1/48频率上加两条数值的虚线。 操作命令如下所示: abline v 1/12, lty "dotted" abline v 1/48, lty "dotted" 执行后输出的结果如下。我们可以看到，SOI数据有一个很明显的主 频率处在1/12上，就是说12个月为一个周期，很明显这是对应总体气温的年度 轮回。然后，我们在下面图上又可以注意到，存在着一个非常不明显的频率。竖 虚线1/48标出了这个频率范围的中心，也就是每四年左右一个循环。前面提到 过，厄尔尼诺现象的周期在3-7年之间。从下面的图中可以看到，厄尔尼诺频率 活动的范围十分地广。厄尔尼诺频率8.运用R语言研究太阳黑子周期频率 下面 我们介绍运用R语言对太阳黑子的周期频率进行估计，从而进一步介绍频率估计 中的频率混叠，即aliasing:。 在我们看电影时，我们常会注意到电影中行驶的车子轮胎向后旋转，这是因为电影胶片的拍摄频率为每秒27张，而真正现实中汽车轮胎的旋转频率远远高于电影胶片的拍摄频率。所以在播放时，汽车轮胎的旋转频率就显得慢了下来，而且方向也被改变了。这种现象叫做频率混叠。 我们这里使用的太阳黑子数据记录是每半年中太阳黑子的数量，即从1749年6月到1978年12月，一共是459个数据，也就是说太阳黑子的数量每一年记录两次。接下来，我们要估计太阳黑子数量的变化周期。一个太阳黑子从生成到消亡可以有几天或者几个月时间。但总体来讲，太阳黑子数量的变化周期大概在11年左右。用上面我们已经介绍过的频率估计程序，我们可以生成下面频率图表。其中主频率在1/22上，因为数据是一年收集两次的，与之对应的周期就是11年一次。下面图表中还显示了一个比较弱的频率在1/140上，这并不是说太阳黑子数量的变化有一个更慢的变化周期，这是因为太阳黑子生成到消亡有几天到几个月的周期，明显地快于数据收集的频率。所以，太阳黑子的生成消亡频率就被混叠了，显示在图表上就变成了一个非常慢的频率。因而在1/140上的频率周期指的是单个太阳黑子生成消亡的周期，而不是太阳黑子数量的变化周期。电脑窗口具体有如下命令: spot scan " spot.per spec.pgram spot, taper 0, log "no" abline v 1/22, lty "dotted" 太阳黑子爆发频率 频率混叠 abline v 1/140, lty "dotted" 参考: 2006, Shumway & Stoffer1