加载中…时间序列分析（Time Series Analysis）
(2012-10-21 16:53:38)
标签：
时间序列rit |
分类： 数据分析 |
时间序列分析（Time Series
Analysis）是指用随机过程理论和数理统计学方法，研究随机数据序列所遵从的统计规律，以用于解决实际问题。由于在多数问题中，随机数据是依时间先后排成序列的，故称为时间序列。时间序列分析在第二次世界大战前就已应用于经济预测。二次大战中和战后，在军事科学、空间科学和工业自动化等部门的应用更加广泛。
一个时间序列通常由4种要素组成：趋势、季节变动、循环波动和不规则波动。
趋势：是时间序列在长时期内呈现出来的持续向上或持续向下的变动。
季节变动：是时间序列在一年内重复出现的周期性波动。它是诸如气候条件、生产条件、节假日或人们的风俗习惯等各种因素影响的结果。
循环波动：是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间，但与趋势不同，它不是朝着单一方向的持续变动，而是涨落相同的交替波动。
不规则波动：是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中，致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列。
时间序列建模基本步骤是：
①用观测、调查、统计、抽样等方法取得被观测系统时间序列动态数据。
②根据动态数据作相关图，进行相关分析，求自相关函数。
③辨识合适的随机模型,进行曲线拟合,即用通用随机模型去拟合时间序列的观测数据。
时间序列分析主要用途：
系统描述：根据对系统进行观测得到的时间序列数据，用曲线拟合方法对系统进行客观的描述。
系统分析：当观测值取自两个以上变量时，可用一个时间序列中的变化去说明另一个时间序列中的变化，从而深入了解给定时间序列产生的机理。
预测未来：一般用ARMA模型拟合时间序列，预测该时间序列未来值。
决策和控制：根据时间序列模型可调整输入变量使系统发展过程保持在目标值上，即预测到过程要偏离目标时便可进行必要的控制。
使用R进行时间序列分析：
一、观察数据
1、获取数据：
>
> jj
[15]
1.000000 0.770000
0.920000 1.000000
1.240000 1.000000
1.160000
[22]
1.300000 1.450000
1.250000 1.260000
1.380000 1.860000
1.560000
[29]
1.530000 1.590000
1.830000 1.860000
1.530000 2.070000
2.340000
[36]
2.250000 2.160000
2.430000 2.700000
2.250000 2.790000
3.420000
[43]
3.690000 3.600000
3.600000 4.320000
4.320000 4.050000
4.860000
[50]
5.040000 5.040000
4.410000 5.580000
5.850000 6.570000
5.310000
[57]
6.030000 6.390000
6.930000 5.850000
6.930000 7.740000
7.830000
[64]
6.120000 7.740000
8.910000 8.280000
6.840000 9.540000
10.260000
[71]
9.540000 8.729999 11.880000
12.060000 12.150000 8.910000 14.040000
[78]
12.960000 14.850000 9.990000 16.200000 14.670000
16.020000 11.610000
2、把数据转变为一个时间序列对象：
> jj = ts(jj,
start=1960, frequency=4)
> jj
1960
0.710000 0.630000
0.850000 0.440000
1961
0.610000 0.690000
0.920000 0.550000
1962
0.720000 0.770000
0.920000 0.600000
1963
0.830000 0.800000
1.000000 0.770000
1964
0.920000 1.000000
1.240000 1.000000
1965
1.160000 1.300000
1.450000 1.250000
1966
1.260000 1.380000
1.860000 1.560000
1967
1.530000 1.590000
1.830000 1.860000
1968
1.530000 2.070000
2.340000 2.250000
1969
2.160000 2.430000
2.700000 2.250000
1970
2.790000 3.420000
3.690000 3.600000
1971
3.600000 4.320000
4.320000 4.050000
1972
4.860000 5.040000
5.040000 4.410000
1973
5.580000 5.850000
6.570000 5.310000
1974
6.030000 6.390000
6.930000 5.850000
1975
6.930000 7.740000
7.830000 6.120000
1976
7.740000 8.910000
8.280000 6.840000
1977
9.540000 10.260000 9.540000
8.729999
1978 11.880000
12.060000 12.150000 8.910000
1979 14.040000
12.960000 14.850000 9.990000
1980 16.200000
14.670000 16.020000 11.610000
这个数据是从1960年开始，个个季度的收入，frequency=4指四个季度。
3、对数据进行绘图：
> plot(jj, ylab="Earnings per Share", main="J &
J")
图 J-J数据的散点图
4、观察数据的季节性变化
> plot(diff(log(jj)),
main="logged and diffed")
图 数据周期性波动图
由上图可以看出，数据随着时间上下波动。
5、使用平滑过滤技巧观察数据的趋势：
在这里，我们用对称的移动平均值来达到过滤和光滑的目的。
下面使用公式：
fjj(t) = ⅛jj(t-2) + ¼jj(t-1) + ¼jj(t) + ¼jj(t+1) +
⅛jj(t+2)
使用如下命令查看趋势图：
> k = c(.5,1,1,1,.5)
> (k = k/sum(k))
> fjj = filter(jj, sides=2, k)
> plot(jj, ylab="Earnings per Share", main="J &
J")
> lines(fjj, col="red")
> lines(lowess(jj), col="blue", lty="dashed")
图 使用平均过滤查看数据趋势图
二、分析数据
1、我们把所有jj数据都取log值。
> log(jj)
2、我们把log值做差，即使用log值数列中第二值减去第一值，第三值减去第二值，第四值减去第三值等等。如果做差处理前数列里有n个数值，处理后的结果中将有n-1个数值。
> dljj = diff(log(jj))
> shapiro.test(dljj)
data:
W = 0.9725, p-value = 0.07211
3、使用柱形分布图和QQ图查看处理后的数据是否符合正态分布
> par(mfrow=c(2,1))
> hist(dljj, prob=TRUE, 12)
> lines(density(dljj))
> qqnorm(dljj)
> qqline(dljj)
http://s3/mw690/83bb57b7tcc8e6085b802&690Series
图 查看数据是否符合正态分布
4、做出延迟图表：
如果现在收集的数据与将来的数据之间存在着正面或者是负面的关联性，我们就可以用现在收集的数据来预测未来的数据。因此找到现在收集到的数据与未来数据之间的关联性是最关键的。找到这种关联性的具体技巧被称作延迟图表，也就是lag.plot。
> lag.plot(dljj, 9, do.lines=FALSE)
图 做出延迟图标
上面语句显示了lag.plot用法，输入的数据是被log和作差后的jj.dat数据。其中9表示我们要考虑从1到9这9个不同的延迟。值得注意的是，在下面图表中显示出延迟4和8显示出了正面关系。其他几个延迟存在着负面关系。
我们可以利用这4和8的正面关系来进行预测，即用现有数据计算接下来的第4个或者是第8个数据。
5、结构拆析。我们把趋势、季节和误差从变型后的jj数据中拆析出来，分别研究，或者分别进行绘图，以便于单独检查。
下面等式代表将要使用的数学模型：
Log(jj) = 趋势 + 季节 + 误差
log(jj) = trend + seasonal + remainder
结构拆析的R命令是stl(), 下面语句中stl命令中输入的是lag变型后的jj数据。其中的“per”输入指的是使用季节循环来进行拆析。stl语句在这里生成了一个叫dog的R物件，然后Plot语句将dog物件进行绘图。
具体操作如下图”
plot(dog
<- stl(log(jj), "per"))
窗口会出现下面所示：
图 结构拆析图
上图中有四行R的绘图，其中第一行代表原来的log（jj）的数据。此数据可以看到总体的上升趋势还存在着一定季节循环性的变化。第二行绘图代表拆析后季节循环的作用。第三行绘图代表拆析后将季节循环作用清除剩余的上升趋势，此数据清楚地看到那种循环性变化已经不存在，剩余的只是趋势。第四行绘图代表将季节循环作用和总体的趋势从数据中清除后所剩余的随机产生的误差。
6、再接下来，我们对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))
> trend = time(jj)-1970
> reg = lm(log(jj)~0+trend+Q, na.action=NULL)
> summary(reg)
Call:
lm(formula = log(jj) ~ 0 + trend + Q, na.action = NULL)
Residuals:
-0.29318 -0.09062 -0.01180 0.08460
0.27644
Coefficients:
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
7、随后，我们将所有估计值插入上面模型中，可得到下面等式：
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")
> lines(fitted(reg), col=2)
输入上面两个命令后，窗口就会显示如下图表。在图表中，黑线代表实际中观察到的log(jj)数据。红线代表由模型计算出来的log(jj)数据。它们看上去非常相似，但区别在于模型中计算出来的数值没误差。
图 实际值与预测值之间的比较图
喜欢
0
赠金笔