• R语言

  • statistics with R 第15章 Time series

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 />
9 天 后
楼主好厉害，辛苦了！ 我顶一个
10 天 后
太好了，请作者继续努力啊，怎么没有后面的内容了呢，
最好最后能够集合在一起，如何能加点自己的心得体会就更爽了，支持！
The traditional spectrum analysis it from Engineering, how every the financial time series are very complicated and have nonlinear term, so other methods are tested, such as fractal Brownian Motion and Hurst statistic for long memory.



However according to my knowledge, the industry does not like too complex ones, and NP is good for VaR, High freq is for fund strategy.



If there is no meaning for arbitrage, no one cares about Hurst nor Fractal BM. Nor Levy, there is a post on Copula, but I think its application is done in market and its econometrics estimation method is too difficult to set up and hardly can find a good theoretical frame work such as ARMA or Kalman filter.
5 天 后
楼主真是好样的，继续努力！

当然最后还要进行总结整理，供大家学习参考！
2 个月 后
1 个月 后
能不能做个文档？？这样方便一点
7 天 后
6 年 后

版主，看到你的翻译很有启发，你可以把其它的翻译发给我吗？还有这本书中的数据从什么地方下载呀？你有的话可以把Statistics with R中的数据发给我吗？

xingzhaoh@163.com，非常感谢