statistics with R 第15章 Time series

本章对比已经考虑过的话题：我们曾对简单统计过程的几种实现十分感兴趣（例如：高斯随机变量、混合高斯过程、线性模型）；现在我们将考虑更复杂过程的实现。

  本章的结构安排如下：

  在引论之后，介绍时间序列的概念并给出几个模拟或者真实的例子，我们将介绍古典的时间序列模型（AR, MA, ARMA, ARIMA, SARIMA），这些模型可以用以刻画时间序列的若干性质。之后介绍时间序列分析的谱方法，用以发现序列中的周期成分。这些模型虽然十分简单具有合理性，但是它们很难完全刻画真实序列的性质：引入基于古典模型的非线性模型（GARCH）。同时引进状态空间模型和kalman滤波：这些模型假设数据建立于线性代数，但我们并不具有所有的信息——有潜在的或者隐含的变量。

其中一些方法很容易可以推广的高维情况，也就是说考虑多变量时间序列，或者说研究同时刻上的几个相关变量，但是此时会出现新的性质（如：协整）。更进一步，如果时间序列的变量太多，则向量模型会出现太多的参数：因此进入面板数据分析领域。

我们最后介绍时间序列的推广：随机过程，此情况下时间为连续的；不规则时间观测序列，此情况下时间是离散的但观测间隔不规则；离散值时间序列（马尔可夫链和隐马尔可夫模型，而不是AR和状态空间模型）。

15.1   引言

15．1．1 例子

  在概率论中，一个时间序列就是（下文你将看到“随机过程”，两者几乎是同样的事情，我们仅会在本章末尾提高随机微分方程时对这两个概念区别对待）一个随机变量序列。在统计学中，时间序列是一个依时间而变的随机变量：例如，股票价格每日变化；空气温度的月观测值；病人每分钟内的平均心率等等。

  plot(LakeHuron,ylab = "",main = "Level of Lake Huron")

  图略

  x <- window(sunspots, start=1750, end=1800)  

  ## window用于截取变量sunspots的1750至1800间的数据（译者注）

  plot(x,ylab = "",main = "Sunspot numbers")

  图略

  你可对上面的序列加上平滑曲线，如：

  plot(x,type ="p",ylab = "",main = "Sunspot numbers")

  k <- 20

  lines( filter(x, rep(1/k,k)),col = "red",lwd = 3 )

  图略

  你可以用竖直柱周期性地下划线一个序列，

data(UKgas)<br />
plot.band <- function (x, ...) {<br />
    plot(x, ...)<br />
    a <- time(x)<br />
    i1 <- floor(min(a))<br />
    i2 <- ceiling(max(a))<br />
    y1 <- par("usr")[3]<br />
    y2 <- par("usr")[4]<br />
    if( par("ylog") ){<br />
        y1 <- 10^y1<br />
        y2 <- 10^y2<br />
    }<br />
    for (i in seq(from=i1, to=i2-1, by=2)){ <br />
        polygon( c(i,i+1,i+1,i),c(y1,y1,y2,y2),col ="grey",border = NA )<br />
    }<br />
    par(new=T)<br />
    plot(x, ...)<br />
    }<br />
plot.band(UKgas,log = "y",ylab = "",main = "UK gas consumption")

 

图略

在时间序列分析之前，我们总是要观察一下这个序列（随时间变化情况，变化类型，混乱值（有时需要检验时间顺序是否正确）；分布性质（直方图，密度估计）；y[i+1]~y[i]，异常值，数据越多异常值可能越多：可能是因为一些人忘记加小数点，也可能是将0或者999代替缺失值等等）。

  x <- LakeHuron

op <- par(mfrow = c(1,2),mar = c(5,4,1,2)+.1,oma = c(0,0,2,0))

hist(x,col = "light blue",xlab = "",main = "")

qqnorm(x,main = "")

qqline(x,col = "red")

par(op)

mtext("Lake Huron levels",line = 2.5,font = 2,cex = 1.2)

  图略

boxplot(x,horizontal = TRUE,col = "pink",main = "Lake Huron levels")

图略

plot(x,ylab = "",main = "Lake Huron levels")

图略

n <- length(x)<br />
k <- 5<br />
m <- matrix(nr=n+k-1, nc=k)<br />
colnames(m) <- c("x[i]", "x[i-1]", "x[i-2]","x[i-3]", "x[i-4]")<br />
for (i in 1:k) { <br />
    m[,i] <- c(rep(NA,i-1), x, rep(NA, k-i))<br />
    }<br />
pairs(m,gap = 0,lower.panel = panel.smooth,<br />
    upper.panel = function (x,y) { panel.smooth(x,y)<br />
    par(usr = c(0, 1, 0, 1))<br />
    a <- cor(x,y, use="pairwise.complete.obs")<br />
    text(.1,.9,<br />
    adj=c(0,1),<br />
    round(a, digits=2),<br />
    col="blue",<br />
    cex=2*a)<br />
    })<br />
title("Lake Huron levels: autocorrelations",line = 3)

图略

我们将会看到其它时间序列的专有图形。第一，自相关函数（ACF）图x[i]与x[i-k]之间的相关系数随k的变化情况（也就上面图形反映的相关关系）。第二，偏自相关函数（PACF）图，它包含同样的信息（详细信息见下文）。最后，谱图，它用以发现序列中的周期成分，频率变化，显示频率变化的重要程度（更多见下文）。

op <- par(mfrow = c(3,1),mar = c(2,4,1,2)+.1)

acf(x, xlab = "")

pacf(x, xlab = "")

spectrum(x, xlab = "", main = "")

par(op)

图略

15．1．2   模拟

  时间序列分析的一个目的（或者步骤）是去探究序列的“结构”，或者说探究序列如何生成，或者说找出可以生成类似序列的算法。下面是一些模拟时间序列。第一个是高斯噪声；第二个是单整噪声（一个随机游走序列）。基于上面两种序列通过加入噪声，线性趋势，周期趋势等等我们可以构造许多其它序列。

<br />
op <- par(mfrow = c(3,3),mar = c(0,0,0,0))<br />
n <- 100<br />
k <- 5<br />
N <- k*n<br />
x <- (1:N)/n<br />
y1 <- rnorm(N)<br />
plot(ts(y1),xlab="", ylab="", main="", axes=F)<br />
box()<br />
y2 <- cumsum(rnorm(N))<br />
plot(ts(y2),xlab="", ylab="", main="", axes=F)<br />
box()<br />
y3 <- cumsum(rnorm(N))+rnorm(N)<br />
plot(ts(y3),xlab="", ylab="", main="", axes=F)<br />
box()<br />
y4 <- cumsum(cumsum(rnorm(N)))<br />
plot(ts(y4),xlab="", ylab="", main="", axes=F)<br />
box()<br />
y5 <- cumsum(cumsum(rnorm(N))+rnorm(N))+rnorm(N)<br />
plot(ts(y5),xlab="", ylab="", main="", axes=F)<br />
box()<br />
# With a trend<br />
y6 <- 1 - x + cumsum(rnorm(N)) + .2 * rnorm(N)<br />
plot(ts(y6),xlab="", ylab="", main="", axes=F)<br />
box()<br />
y7 <- 1 - x - .2*x^2 + cumsum(rnorm(N)) +.2 * rnorm(N)<br />
plot(ts(y7),xlab="", ylab="", main="", axes=F)<br />
box()<br />
# With a seasonnal component<br />
y8 <- .3 + .5*cos(2*pi*x) - 1.2*sin(2*pi*x) +.6*cos(2*2*pi*x) + .2*sin(2*2*pi*x) +-.5*cos(3*2*pi*x) + .8*sin(3*2*pi*x)<br />
plot(ts(y8+ .2*rnorm(N)),xlab="", ylab="", main="", axes=F)<br />
box()<br />
lines(y8, type="l", lty=3, lwd=3, col="red")<br />
y9 <- y8 + cumsum(rnorm(N)) + .2*rnorm(N)<br />
plot(ts(y9),<br />
xlab="", ylab="", main="", axes=F)<br />
box()<br />
par(op)<br />

图略

15．1．3    时间序列分析的主要问题

  在统计学中，数据的独立性是我们希望的性质，问题是时间序列数据间往往是相关的。时间序列分析的目的就是要通过统计模型提取一个序列的结构，将原始序列转化为相互独立序列。我们可以从另一个角度来看这个问题：当你研究一个统计现象时，你往往有多个观测值，而当你研究时间序列问题时，你仅有一次观测值。因此我们可以用不同时点一次观测值来代替同一时点上的多个观测值。依统计问题不同，这两点是否能够等同起来——这个问题也就是所谓的遍历性——更多见下文。

15．1．4   自相关

  由于时间序列观测值往往不独立，因此我们首先可以看看一下它们的相关性：滞后k期的自相关（ACF）系数为（严格来说我们可以从样本观测值中考察样本ACF，或者我们由理论模型算得理论自相关函数）是n期与n-k期观测值之间的相关系数。为了由样本中算得ACF，必须假设它不依赖于样本观测期n，仅依赖于滞后期k。我们可以由下面得函数计算ACF。

<br />
my.acf <- function (x,lag.max = ceiling(5*log(length(x)))) {<br />
m <- matrix(c( NA,rep( c(rep(NA, lag.max-1), x),lag.max ),rep(NA,, lag.max-1)),byrow=T,nr=lag.max)<br />
    x0 <- m[1,]<br />
    apply(m,1,cor, x0, use="complete")<br />
    }<br />
n <- 200<br />
x <- rnorm(n)<br />
plot(my.acf(x),xlab = "Lag",type ="h")<br />
abline(h=0)<br />

图略

事实上，系统已经有一个这样的函数，

<br />
op <- par(mfrow=c(2,1))<br />
acf(x, main="ACF of white noise")<br />
x <- LakeHuron<br />
acf(x, main="ACF of a time series (Lake Huron)")<br />
par(op)<br />
 

图略

这是引言中模拟数值的自相关函数，该图形是托尾的：你能从中判别正确的滞后阶数吗？

##原书此处代码可能有误

依我来看，可以看到三组：当ACF几乎为零时，数据不相关；当ACF时正时负时，数据也许是周期的；另一种情况，数据是相关的。

15．1．5   白噪声

一个白噪声序列是一个不相关随机变量序列，它们的期望为零，方差为常数。换句话说，它们是独立同分布随机变量，至二阶矩（它们也许是相关的，但为不能用线性相关来考察的非线性相关关系；它们也许是不同分布的，只要它们的均值和方差相同）。我们经常费尽心机要把我们的序列分解成一个“趋势项”加噪声项，这个噪声项也许是白噪声的。例如一个独立同分布随机变量序列是一个白噪声序列。

<br />
my.plot.ts <- function (x, main="") { <br />
    op <- par(mar=c(2,2,4,2)+.1)<br />
    layout( matrix(c(1,2),nr=1,nc=2), widths=c(3,1) )<br />
    plot(x, xlab="", ylab="")<br />
    abline(h=0, lty=3)<br />
    title(main=main)<br />
    hist(x, col="light blue", main="", ylab="", xlab="")<br />
    par(op)<br />
    }<br />
n <- 100<br />
x <- ts(rnorm(n))<br />
my.plot.ts(x, "Gaussian iid noise")<br />

图略

一个独立同分布零均值随机变量序列也是一个白噪声序列。

<br />
n <- 100<br />
x <- ts(runif(n,-1,1))<br />
my.plot.ts(x, "Non gaussian iid noise")<br />

图略

<br />
x <- ts(rnorm(100)^3)<br />
my.plot.ts(x, "Non gaussian iid noise")<br />

图略

也有不独立但不相关序列：根据定义不相关仅为不独立至二阶矩。

<br />
n <- 100<br />
x <- rep(0,n)<br />
z <- rnorm(n)<br />
for (i in 2:n) { <br />
    x[i] <- z[i] * sqrt( 1 + .5 * x[i-1]^2 )<br />
    }<br />
my.plot.ts(x, "Non iid noise")<br />

图略

一些确定序列也看起来象白噪声序列

<br />
n <- 100<br />
x <- rep(.7, n)<br />
for (i in 2:n) { <br />
    x[i] <- 4 * x[i-1] * ( 1 - x[i-1] )<br />
}<br />
op <- par()<br />
layout( matrix(c(1,2),nr=1,nc=2), widths=c(3,1) )<br />
x <- x - mean(x)<br />
plot(ts(x))<br />
abline(h=0, lty=3)<br />
title(main="A deterministic time series")<br />
hist(x, col="light blue", main="", ylab="")<br />
par(op)<br />

图略

<br />
n <- 1000<br />
tn <- cumsum(rexp(n))<br />
# A C^infinity function defined as a sum<br />
# of gaussian densities<br />
f <- function (x) { <br />
    # If x is a single number: sum(dnorm(x-tn))<br />
apply( dnorm( outer(x,rep(1,length(tn))) –<br />
outer(rep(1,length(x)),tn) ),1,sum )<br />
    }<br />
op <- par(mfrow=c(2,1))<br />
curve(f(x),xlim = c(1,500),n = 1000,main = "From far away, it looks random...")<br />
curve(f(x),xlim = c(1,10),n = 1000,main="...but it is not: it is a C^infinity function")<br />
par(op)<br />

我们有一些检验方法可以检验一个序列是否为白噪声。

15．1．6   诊断：是白噪声序列吗？

时间序列分析主要是寻找一种方法来拟和这个序列（或者就像我们引言所讲，由白噪声序列来生成一个近似的序列）。我们可以以如下方向来建立模型：考察数据序列并将其转换成类似白噪声序列。为了检验我们的分析是否正确，需要检验模型残差是否为白噪声序列。至此，我们可以考察自相关函数图（平均而言，仅能有5％的相关系数线超过虚线，如果有更多，那么我们的分析或者说结果是有疑问的）。

<br />
z <- rnorm(200)<br />
op <- par(mfrow=c(2,1), mar=c(5,4,2,2)+.1)<br />
plot(ts(z))<br />
acf(z, main = "")<br />
par(op)<br />

图略

<br />
x <- diff(co2)<br />
y <- diff(x,lag=12)<br />
op <- par(mfrow=c(2,1), mar=c(5,4,2,2)+.1)<br />
plot(ts(y))<br />
acf(y, main="")<br />
par(op)<br />

图略



  为了得到数值结果（p值），我们可以进行Box－Pierce或者Ljung－Box 检验（这些检验也称为“portmanteau 检验”）：其思想是考察（加权）平均和一阶自相关系数，这些求和结果（渐近）服从卡方分布（Ljung－Box 检验较Box－Pierce检验更为稳健，小样本性质更好近似的卡方分布）。

> Box.test(z) # Box-Pierce

Box-Pierce test

X-squared = 0.014, df = 1, p-value = 0.9059

> Box.test(z, type="Ljung-Box")

Box-Ljung test

X-squared = 0.0142, df = 1, p-value = 0.9051

> Box.test(y)

Box-Pierce test

X-squared = 41.5007, df = 1, p-value = 1.178e-10

> Box.test(y, type=”Ljung”)

Box-Ljung test

X-squared = 41.7749, df = 1, p-value = 1.024e-10

<br />
op <- par(mfrow=c(2,1))<br />
plot.box.ljung <- function (z,k = 15,main = "p-value of the Ljung-Box test",ylab = "p-value") {<br />
    p <- rep(NA, k)<br />
    for (i in 1:k) { <br />
        p[i] <- Box.test(z, i,<br />
        type = "Ljung-Box")$p.value<br />
    }<br />
    plot(p,type = "h",ylim = c(0,1),lwd = 3,main = main,ylab = ylab)<br />
    abline(h = c(0,.05),lty = 3)<br />
    }<br />
plot.box.ljung(z, main="Random data")<br />
plot.box.ljung(y, main="diff(diff(co2),lag=12)")<br />
par(op)<br />

图略

还有其它的检验方法：McLeod-Li, Turning-point, difference-sign, rank 检验等等。我们也可以用DW检验，在回归问题中我们已经应用了DW检验。

TODO: check that I actually mention it.

library(car)

?durbin.watson

library(lmtest)

?dwtest

这里应用的同样的检验方法，但并非以同样的方式检验，结果也许不同：

> dwtest(LakeHuron ~ 1)

Durbin-Watson test

data: LakeHuron ~ 1

DW = 0.3195, p-value = < 2.2e-16

alternative hypothesis: true autocorrelation is greater than 0

> durbin.watson(lm(LakeHuron ~ 1))

lag Autocorrelation D-W Statistic p-value

1 0.8319112 0.3195269 0

Alternative hypothesis: rho != 0

<br />
op <- par(mfrow=c(2,1))<br />
library(lmtest)<br />
plot(LakeHuron,main = "Lake Huron")<br />
acf(LakeHuron,main = paste("Durbin-Watson: p =",signif( dwtest( LakeHuron ~ 1 ) $ p.value, 3 )))<br />
par(op)<br />

图略

<br />
n <- 200<br />
x <- rnorm(n)<br />
op <- par(mfrow=c(2,1))<br />
x <- ts(x)<br />
plot(x, main="White noise", ylab="")<br />
acf(x,main = paste("Durbin-Watson: p =",signif( dwtest( x ~ 1 ) $ p.value, 3)))<br />
par(op)<br />

图略

<br />
n <- 200<br />
x <- rnorm(n)<br />
op <- par(mfrow=c(2,1))<br />
y <- filter(x,.8,method="recursive")<br />
plot(y, main="AR(1)", ylab="")<br />
acf(y,main = paste("p =",signif( dwtest( y ~ 1 ) $ p.value, 3 )))<br />
par(op)<br />

图略

但是注意：默认的检验是检验自相关是零还是正的，如果自相关显著为负，结果会令人费解。

<br />
set.seed(1)<br />
n <- 200<br />
x <- rnorm(n)<br />
y <- filter(x, c(0,1), method="recursive")<br />
op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)<br />
plot(y,main = paste("one-sided DW test: p =",signif( dwtest ( y ~ 1 ) $ p.value, 3 )))<br />
acf( y, main="")<br />
pacf(y, main="")<br />
par(op)<br />

图略

<br />
op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)<br />
res <- dwtest( y ~ 1, alternative="two.sided")<br />
plot(y,main = paste("two-sided p =",signif( res$p.value, 3 )))<br />
acf(y, main="")<br />
pacf(y, main="")<br />
par(op)<br />

图略

还有其它检验，例如runs检验（它看起来像数据在run，一个run是对同一个符号的连续观测；这类检验主要用来进行时间序列的定性检验）。

library(tseries)

?runs.test

或者Cowles-Jones 检验

TODO

A sequence is a pair of consecutive returns of the same sign; a

reversal is a pair of consecutive returns of opposite signs.

Their ratio, the Cowles-Jones ratio,

number of sequences

CJ = ---------------------

number of reversals

should be around one IF the drift is zero.

15．1．7    tsdiag函数

  事实上，系统已经有一个函数，“tsdiag”，这个函数可以执行若干时间序列的检验（plot， ACF， Ljung－Box检验）。

<br />
data(co2)<br />
r <- arima(co2,order = c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12))<br />
tsdiag(r)<br />

图略

15．2    简单时间序列模型

15．2．1    古典模型

   一般地，我们常费尽心机将一个时间序列分解成三个部分：一个趋势成分（常常是一个仿真函数），一个季节成分（周期函数）和白噪声项。本节中，我们将介绍经典地时间序列模型（主要是回归），其中一些方法是相互关联的、有用的，但是也有一些仅仅是思想而已，实用性不大。上眼吧！

本节的结构安排如下：我们首先用回归方法建立一个时间序列模型，具体的是一个多项式函数加上一个正弦波，或者更一般地，一个多项式函数加上一个周期信号。然后我们介绍基于（指数）移动平均地建模技术，并假设序列看起来像一个常数项加上一个噪声，或者一个线性函数加上一个噪声项，也就是所谓地HW滤波。

15．2．2   第一次尝试：回归

<br />
data(co2)<br />
plot(co2)<br />

图略

此处，我们假设回归时，时间项为解释变量，现在让我们将这组数据写成是多项式函

数与正玹函数地和。

<br />
y <- as.vector(co2)<br />
x <- as.vector(time(co2))<br />
r <- lm( y ~ poly(x,1) + cos(2*pi*x) + sin(2*pi*x) )<br />
plot(y~x, type="l", xlab="time", ylab="co2")<br />
lines(predict(r)~x, lty=3, col="red", lwd=3)<br />

图略

趋势项仅仅是一次多项式也许是不够地，虽然观察预测点还不能清楚的看到这一点，

但是让我们看一下残差，结论就一幕了然了。

<br />
plot( y-predict(r),main = "The residuals are not random yet",xlab = "Time",ylab = "Residuals" )<br />

图略

让我把模型变得更加复杂些：二次多项式加正玹函数（如果二次多项式还不够，就就

用样条曲线代替多项式）。

<br />
r <- lm( y ~ poly(x,2) + cos(2*pi*x) + sin(2*pi*x) )<br />
plot(y~x, type="l", xlab="Time", ylab="co2")<br />
lines(predict(r)~x, lty=3, col="red", lwd=3)<br />

图略

<br />
plot( y-predict(r),main = "Better residuals -- but still not random",xlab = "Time",ylab = "Residuals" )<br />

图略

我们也可以重新修正周期成分，但是可能这并不会改善模型效果。

<br />
r <- lm( y ~ poly(x,2) + cos(2*pi*x) + sin(2*pi*x)+ cos(4*pi*x) + sin(4*pi*x) )<br />
plot(y~x, type="l", xlab="Time", ylab="co2")<br />
lines(predict(r)~x, lty=3, col="red", lwd=3)<br />

图略

<br />
plot( y-predict(r),type = "l",xlab = "Time",ylab = "Residuals",main = "Are those residuals any better?" )<br />

图略

然而，通过ACF图我们可以看到，这些模型并不那么如我们所期望的那样好：如果残

差真的是白噪声的，你们ACF图中，将很少有自相关线超过虚线…

<br />
r1 <- lm( y ~ poly(x,2) +cos(2*pi*x) +sin(2*pi*x) )<br />
r2 <- lm( y ~ poly(x,2) +cos(2*pi*x) +sin(2*pi*x) +cos(4*pi*x) +sin(4*pi*x) )<br />
op <- par(mfrow=c(2,1))<br />
acf(y - predict(r1))<br />
acf(y - predict(r2))<br />
par(op)<br />

图略

由这些图形可以得到以下两点事实：第一，毕竟修正周期成分是有益的；第二，模型

仍然存在残差自相关的问题。

15.2.3   其它尝试（明显错误的想法）

为了估计周期成分，我们可以取一月份观测值的平均数，二月份观测值的平均数等

等。

<br />
m <- tapply(co2, gl(12,1,length(co2)), mean)<br />
m <- rep(m, ceiling(length(co2)/12)) [1:length(co2)]<br />
m <- ts(m, start=start(co2), frequency=frequency(co2))<br />
op <- par(mfrow=c(3,1), mar=c(2,4,2,2))<br />
plot(co2)<br />
plot(m, ylab = "Periodic component")<br />
plot(co2-m, ylab = "Withou the periodic component")<br />
r <- lm(co2-m ~ poly(as.vector(time(m)),2))<br />
lines(predict(r) ~ as.vector(time(m)), col=¡¯red¡¯)<br />
par(op)<br />

图略

然而，观察一下这个模型的残差就会发现，周期成分依然存在，

<br />
op <- par(mfrow=c(4,1), mar=c(2,4,2,2)+.1)<br />
plot(r$res, type = "l")<br />
acf(r$res, main="")<br />
pacf(r$res, main="")<br />
spectrum(r$res, col=par("fg"), main="")<br />
abline(v=1:6, lty=3)<br />
par(op)<br />

图略

15.2.4    其它尝试（效果明显比以前要好）

用同期均值估计量估计趋势成分仅当采用如下形式的模型时是有效的，

y ~ a + b x + periodic component

但是这里采用二次项，你们模型就无效了。然而我们可以重新采用同样的简单回归方法找寻周期成分和趋势成分。这里周期成分可以是任意的：在上面的模型中有12项系数需要估计。

这里采用的方法是虚拟变量方法，在用二次多项式拟和趋势成分的同时，用12个虚拟

变量反映周期变化因素（译者注）。

<br />
k <- 12<br />
m <- matrix( as.vector(diag(k)),nr = length(co2),nc = k,byrow = TRUE )<br />
m <- cbind(m, poly(as.vector(time(co2)),2))<br />
r <- lm(co2~m-1)<br />
summary(r)<br />
b <- r$coef<br />
y1 <- m[,1:k] %*% b[1:k]<br />
y1 <- ts(y1,start=start(co2),frequency=frequency(co2))<br />
y2 <- m[,k+1:2] %*% b[k+1:2]<br />
y2 <- ts(y2,start=start(co2),frequency=frequency(co2))<br />
res <- co2 - y1 - y2<br />
op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)<br />
plot(co2)<br />
lines(y2+mean(b[1:k]) ~ as.vector(time(co2)),col="red")<br />
plot(y1)<br />
plot(res)<br />
par(op)<br />

图略

这里的分析并没有完成，我们仍然需要研究一下模型残差——从中可以看出我们已经去掉的周期成分。

<br />
op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)<br />
acf(res, main="")<br />
pacf(res, main="")<br />
spectrum(res, col=par("fg"), main="")<br />
abline(v=1:10, lty=3)<br />
par(op)<br />

图略

15．    2.5    一些思想，样条曲线

在上一个例子中，如果周期更长，我们可以用样条曲线寻找平滑周期成分，并且只有更

少的估计参数：最多15个。

TODO

Exercise left to the reader...

15．2．6    寻找趋势（或者去除季节成分）：移动平均

要寻找一个序列的趋势成分，你平滑它就可以了，例如移动平均平滑（MA）。

<br />
x <- co2<br />
n <- length(x)<br />
k <- 12<br />
m <- matrix( c(x, rep(NA,k)), nr=n+k, nc=k )<br />
y <- apply(m, 1, mean, na.rm=T)<br />
y <- y[1:n + round(k/2)]<br />
y <- ts(y, start=start(x), frequency=frequency(x))<br />
y <- y[round(k/4):(n-round(k/4))]<br />
yt <- time(x)[ round(k/4):(n-round(k/4)) ]<br />
plot(x, ylab="co2")<br />
lines(y~yt, col="red")    <br />

图略

移动平均平滑得到了确定的函数项成分（例如，用直到三阶的多项式——选择合适的系数）

系统“filter”函数也具有上面的功能。

<br />
x <- co2<br />
plot(x, ylab="co2")<br />
k <- 12<br />
lines( filter(x, rep(1/k,k)), col="red")<br />

最近比较忙，我会把本章翻译完