时间序列分析
需要解释的概念

平稳时间序列

平稳时间序列
平稳时间序列及其数字特征


均值为常数，协方差函数与时间无关的时间序列叫做平稳时间序列


平稳时间序列有一种简化模型，就是均值为0的平稳时间序列。由于平稳时间序列的均值是常数，所以当某个平稳时间序列的均值非0的时候，可以给每一个元素都减去他们的均值，而且这个均值是指每个随机变量的均值，而不是指样本序列的均值，所以只要知道随机变量是什么可以很容易得到他的均值，因此任意一个平稳时间序列都可以转换为零均值的平稳时间序列。


均值讨论了之后，本来应该讨论自相关函数，但是不知道为什么，这里重要的数字特征引入了一个自相关系数函数（Autocorrelatrion Coefficient Function，ACF），就是那个皮尔逊相关系数的那个相关系数，计算也是协方差除以两个随机变量的方差之积开根号。非常神奇。这里引入一个计算公式：rou_k = v_k/v_0，这里rou_k就是自相关系数函数ACF，v_k是X_t和X_t+k的协方差函数，v_0同理，但是因为k=0我们倾向于叫v_0为方差。


平稳时间序列最重要的两个数字特征应该就是均值，协方差和自相关系数函数ACF（名字有点绕口我得多打几遍才能记住），由于均值比较简单，而且一般我们研究的都是零均值的平稳时间序列，因此我们下面讨论协方差函数v_k和自相关系数函数ACF的性质。


v_k和rou_k都是偶函数。这也比较好理解，因为根据他们的定义，k取一个相反数实际上也就相当于求协方差的时候或者求相关系数的时候两个参数调换一下位置。


v_k的模始终小于v_0，rou_k的模始终小于1.ACF始终小于1可以理解，毕竟相关系数嘛。因为协方差其实跟相关系数差不多，反应的其实也都是两个随机变量之间的相关关系，两个随机变量之间的相关关系再相关也不会比自己和自己更相关了吧，所以协方差函数大小始终小于方差函数大小。


非负定性。其实一直我也不老理解非负定性的概念的，大概查了查资料，大概懂了是说某个矩阵的ij位置元素作为前乘向量的第i个元素和后乘向量第j个元素乘积的系数的二次多项式，如果这个二次多项式是非负的，那么这个矩阵就具有非负定性。其实无论是v_k还是rou_k都应该是一个矩阵的，但是因为他们都跟时间无关，因此就自动认为把他们跟时间没关系了，可以用向量来表示。好像是因为把z_i跟X_t-i点积了一下然后求方差非负，可以化成协方差矩阵的非负定性。然后因为协方差v_k除以v_0直接就是rou_k所以也就得到rou_k对应的矩阵具有非负定性。


以上数字特征就介绍的差不多了，整个平稳时间序列也介绍的差不多了。理论上的东西讲完了之后，要看一下样本里面怎么处理。


样本均值X_bar是样本观测值的平均值，样本协方差函数v_k_hat是1/n乘以t=1到t=(n-k)求和，(X_t - X_bar)(X_t+k - X_bar)，样本自相关系数协方差函数rou_k_hat是v_k_hat除以v_0_hat。


这里面样本均值和样本自相关系数函数都没啥可说的，就是正常情况，这里要提一下样本协方差函数。有些人可能会问，明明求和只有n-k项，为什么要除以1/n？这个问题好像涉及到一点正定性的原因。我找到了一份参考材料，在这份参考材料中提到，采取1/n的劣势是如果k过大的话，整体误差会增大，但是由于我们一般不用太大的k计算v_k_hat，而且更重要的是，只有1/n的估计式才是正定的。具体情况可以再仔细看看参考资料，这部分也不是很重要，就不再仔细解释了。


样本的部分基本上也就介绍这些了，接下来是一种特殊的平稳时间序列：白噪声序列


零均值，协方差函数满足方差为sigma^2以外其余都等于0的平稳时间序列就是白噪声序列。如果每个随机变量都服从正态分布的话，这个白噪声序列被称为正态白噪声序列。


由于在实际上，正态白噪声序列比较难以判断，因此引入一个定理，可以大致根据序列的相关函数对是否是正态白噪声序列进行初步判断。


对于正态白噪声序列来说，一切k>0，rou_k_hat服从N(0, 1/n)。因此，根据这个性质，rou_k_hat落在两倍标准差内的概率为95.5%，可以根据这个性质大致判断。


上面部分，最基本的概念就引入完毕了，接下来要进入重头戏了。时间序列分析的基本内容之一就是通过分析ACFrou_k等数字特征来研究X_t和X_t的滞后项（即X_t-1,X_t-2等）之间的线性相关关系，并且选取适当的模型来刻画这种关系，从而可以用X_t的样本观测值来对时间序列在未来某些时刻的值进行预报等。为了描述方便，此时引入延迟算子的概念。


BX_t = BX_t-1，则称B为一步延迟算子，B^k称为k步延迟算子。


由于上面说到了要研究X_t与滞后项的关系，因此引入一个随机变量序列Y_t = 求和i=0到p a_i*X_t-i，注意这里面有一个p作为参数，也就是向前追溯多少步。这个随机变量序列也被称为X_t的滑动加权平均序列。


介绍了延迟算子自然要用上，Y_t可以用延迟算子的函数和X_t来表示。即Y_t = f(B)X_t。f(B)是B的p次多项式，为 求和i=0到p a_iB^i。


大概总结一下，如下图所示

平稳时间序列分析
方法性工具
差分运算

一阶差分
p阶差分
k步差分

延迟算子
线性差分方程
ARMA模型
AR模型


AR模型平稳性判别

特征根判别
平稳域判别



平稳AR模型的统计性质


MA模型


MA模型的统计性质

均值
方差
自协方差q阶截尾
自相关系数q阶截尾



MA模型的可逆性


偏自行管系数截尾


ARMA模型

平稳条件 & 可逆条件
传递形式 & 逆转形式
ARMA(p,q)模型的统计性质

平稳序列建模
建模步骤

求序列的样本自相关系数ACF和样本偏自相关系数PACF值
根据两个系数的性质，选择阶数适当的ARMA(p,q)模型进行拟合
估计模型中未知参数的值
检验模型有效性。如果拟合不符，返回步骤2重新选择模型拟合
模型优化。通过检验，返回步骤2，充分考虑各种可能，建立多个拟合模型，选择最优模型
利用拟合模型，预测序列将来走势

样本自相关系数 & 偏自相关系数
模型识别
参数估计

矩估计
极大似然估计
最小二乘估计

模型检验

模型显著性检验
参数显著性检验

模型优化

问题提出
AIC准则
SBC（BIC）准则

序列预测
线性预测函数
预测方差最小原则
线性最小方差与预测的性质

条件无偏最小方差估计值
AR§序列预测
MA(q)序列预测
ARMA(p,q)序列预测


#AR模型平稳性判别
#1、arima.sim函数拟合
#arima.sim(n,list(ar=,ma=,order=),sd=)
#n：拟合序列长度
#list：指定具体的模型参数
#（1）拟合平稳AR(p)模型，要给出自回归系数，如果指定拟合的AR模型为非平稳模型，系统会报错
#（2）拟合MA(q)模型，要给出移动平均系数
#（3）拟合平稳ARMA(p,q)模型需要同时给出自回归系数和移动平均系数，如果指定拟合的ARMA模型为非平稳模型，系统会报错
#（4）拟合ARIMA(p,d,q)模型（第5章介绍），除了需要给出自回归系数和移动平均系数，还需要增加order选项。order=c(p,d,q)
# 其中，p为自回归阶数；d为差分阶数；q为移动平均阶数
#sd：指定序列的标准差，不特殊指定，系统默认 sd=1

#2、filter函数拟合
#filter(e,filter=,method=,circular=)
#e：随机波动序列的变量名
#filter：指定模型系数
#（1）AR(p)模型为filter=c(∅1,∅2,…,∅p)
#（2）MA(q)模型为filter=c(1,-θ1,-θ2,…,-θq)
#method：指定拟合的是AR模型还是MA模型
#（1）method=“recursive”为AR模型
#（2）method=“convolution”为MA模型
#circular：拟合MA模型时专用的一个选项，circular=T可以避免NA数据出现

#3.1
x1=arima.sim(n=100,list(ar=0.8))
x3=arima.sim(n=100,list(ar=c(1,-0.5)))
e=rnorm(100)
x2=filter(e,filter=-1.1,method = "recursive")
x4=filter(e,filter = c(1,0.5),method="recursive")

layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))

#layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
# layout() 输入一个矩阵，2行2列，
#第1个图占两个1，第2个图占2，第3个图占3

ts.plot(x1) #平稳
ts.plot(x2) #非平稳
ts.plot(x3) #平稳
ts.plot(x4) #非平稳

#3.5
x1=arima.sim(n=1000,list(ar=0.8))
x2=arima.sim(n=1000,list(ar=-0.8))
x3=arima.sim(n=1000,list(ar=c(1,-0.5)))
x4=arima.sim(n=1000,list(ar=c(-1,-0.5)))

#AR模型样本自相关图
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
acf(x1)
acf(x2)
acf(x3)
acf(x4)

#AR模型样本偏自相关图
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)

#3.6
x1=arima.sim(n=1000,list(ma=-2))
x2=arima.sim(n=1000,list(ma=-0.5))
x3=arima.sim(n=1000,list(ma=c(-4/5,16/25)))
x4=arima.sim(n=1000,list(ma=c(-5/4,25/16)))

layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
acf(x1)
acf(x2)
acf(x3)
acf(x4)

layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)

#3.8
x=arima.sim(n=1000,list(ar=0.5,ma=-0.8))
acf(x)
pacf(x)

#3.9
#读入数据，并绘制时序图
layout(matrix(c(1), 1, byrow = TRUE))

a=read.table("E:/data/file8.csv",sep=",",header=T)
x=ts(a$kilometer,start=1950)
plot(x)

#白噪声检验
for(i in 1:2) print(Box.test(x,type="Ljung-Box",lag=6*i))

#绘制自相关图和偏自相关图
acf(x)
pacf(x)

#3.10（有误）
#读入数据，并绘制时序图
overshort<-read.table("E:/data/file9.csv",sep=",",header=T)
overshort<-ts(overshort)
plot(overshort)

#白噪声检验
for(i in 1:2) print(Box.test(overshort,type="Ljung-Box",lag=6*i))

#绘制自相关图和偏自相关图
acf(overshort)  #一阶截尾
pacf(overshort)  #拖尾

#3.11
#读入数据，并绘制时序图
b=read.table("E:/data/file10.csv",sep=",",header=T)
dif_x=ts(diff(b$change_temp),start=1880)
plot(dif_x)

#白噪声检验
for(i in 1:2) print(Box.test(dif_x,type="Ljung-Box",lag=6*1))

#绘制自相关图和偏自相关图
acf(dif_x)
pacf(dif_x)

#auto.arima函数：先安装packages:zoo & forecast，然后用library调用程序包
#auto.arima(x,max.p=5,max.q=,ic=)
#x：需要定阶的序列名
#max.p：自相关系数最高阶数，不特殊指定的话，系统默认值为5
#max.q：移动平均系数最高阶数，不特殊指定的话，系统默认值为5
#ic：指定信息量准则，ic有"aicc","aic"和"bic"三个选项，系统默认AIC准则

library(zoo)
library(forecast)

#3.9 系统自动定阶
auto.arima(x)

#3.10 系统自动定阶
auto.arima(dif_x)

#3.11 系统自动定阶
auto.arima(dif_x)

#arima(x,order=,include.mean=,method=)
#x：要进行模型拟合的序列名
#order：指定模型阶数。order=c(p,d,q)
#（1）p为自回归阶数
#（2）d为差分阶数，本章不涉及缠粉问题，所以d=0
#（3）q为移动平均阶数
#include.mean：要不要包含常数项
#（1）include.mean=T，需要拟合常数项，这也是系统默认设置
#（2）如果不需要拟合常数项，需要特别指定include.mean=F
#method：指定参数估计方法
#（1）method="CSS-ML"，默认的是条件最小二乘与极大似然估计混合方法
#（2）method="ML"，极大似然估计
#（3）mrthod="CSS"，条件最小二乘估计

#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header=T)
x=ts(a$kilometer,start=1950)
x.fit=arima(x,order=c(2,0,0),method="ML")
x.fit

#3.10续（有误）
overshort=read.table("E:/data/file9.csv",sep=",",header = T)
overshort=ts(overshort)
overshort.fit=arima(overshort,order=c(0,0,1))
overshort.fit

#3.11续
b=read.table("E:/data/file10.csv",sep=",",header = T)
dif_x=ts(diff(b$change_temp),start = 1880)
dif_x.fit=arima(dif_x,order=c(1,0,1))
dif_x.fit

#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0),method = "ML")
for(i in 1:2) print(Box.test(x.fit$residuals,lag=6*i))

#3.10续（有误）
overshort=read.table("E:/data/file9.csv",sep=",",header = T)
overshort=ts(overshort)
overshort.fit=arima(overshort,order=c(0,0,1))
for(i in 1:2) print(Box.test(overshort.fit$residual,lag=6*i))

#3.11续
b=read.table("E:/data/file10.csv",sep=",",header = T)
dif_x=ts(diff(b$change_temp),start = 1880)
dif_x.fit=arima(dif_x,order=c(1,0,1),method="CSS")
for(i in 1:2) print(Box.test(dif_x.fit$residual,lag=6*i))

#pt(t,df=,lower.tail=)
#t：t统计量的值
#df：自由度
#lower.tail：确定计算概率的方向
#（1）lower.tail=T，计算Pr(X≤x)。对于参数显著性检验，如果参数估计值为负，选择lower.tail=T
#（2）lower.tail=F，计算Pr(X>x)。对于参数显著性检验，如果参数估计值为正，选择lower.tail=F

#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0),method="ML")
x.fit

#3.15（有误）
#读入数据，绘制时序图
x=read.table(file="E:/data/file11.csv",sep=",",header = T)
x=ts(x)
plot(x)

#序列白噪声检验
for(i in 1:2) print(Box.test(x,lag=6*i))

#绘制自相关图和偏自相关图
acf(x)
pacf(x)

#拟合MA(2)模型
x.fit1=arima(x,order=c(0,0,2))
x.fit1

#MA(2)模型显著性检验
for(i in 1:2) print(Box.test(x,fit1$residual,lag=6*1))

#拟合AR(1)模型
x.fit2=arima(x,order=c(1,0,0))
x.fit2

#AR(1)模型显著性检验
for(i in 1:2) print(Box.test(x.fit2$residual,lag=6*i))


#forecast(objuct,h=,level=)
#object：拟合信息文件名
#h：预测期数
#level：置信区间的置信水平。不特殊指定的话，系统会自动给出置信水平分别为80%和95%的双层置信区间

#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0))
x.fore=forecast(x.fit,h=5)
x.fore

#系统默认输出预测图
plot(x.fore)

#个性化输出预测图
L1=x.fore$fitted-1.96*sqrt(x.fit$sigma2)
U1=x.fore$fitted+1.96*sqrt(x.fit$sigma2)
L2=ts(x.fore$lower[,2],start = 2009)
U2=ts(x.fore$upper[,2],start = 2009)
c1=min(x,L1,L2)
c2=max(x,L2,U2)
plot(x,type="p",pch=8,xlim=c(1950,2013),ylim=c(c1,c2))
lines(x.fore$fitted,col=2,lwd=2)
lines(x.fore$mean,col=2,lwd=2)
lines(L1,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U2,col=4,lty=2)

文章目录
一、时间序列与平稳序列
1.时间序列的概念
2.重要的时间序列——平稳序列
3.特殊的平稳序列——白噪声
4.多平稳序列的相互关系
回顾总结

一、时间序列与平稳序列
1.时间序列的概念
时间序列，就是按照时间次序排列的随机变量列，其最重要的特征就是具有时间关系，即处于不同时间的随机变量可能具有一定的联系。生活中有许多时间序列，如每个月的平均气温、股市每天的收盘价等等，都是时间序列。
任何时间序列，经过合理的变换后都可以看作由三个部分叠加而成：趋势项，周期项和随机噪声项。趋势项大体刻画了时间序列的变化趋势，是一个固定的、可以预测的项；周期项是具有一定周期的时间序列，比如一年四季每个季节有各自的特征，就可以用周期项来刻画；随机噪声则是随机干扰，一般被视为独立的零均值序列。
以上关系概括说来，就是

$X_t=T_t+S_t+R_t.$

在实际生活中，时间不能倒流，所以时间序列往往只能够发生一次，即获得一次观测。$X_1,X_2,\cdots$的一组实际数值$x_1,x_2,\cdots$是时间序列的一次实现或一条轨道。
在获得观测值后，要对时间序列进行以上的分解，才能够获得具有实际意义的分布。有一些常用的分解方式，如分段趋势分解，回归直线法，二次曲线回归法，逐步平均法等等。
随机过程中将时间指标分成连续集与离散集两种，即$t$的取值可以是连续的$\R,\R_+$或离散的$\Z,\Z_+$，我们将重点放在离散时间序列上。
2.重要的时间序列——平稳序列
时间序列的趋势项和季节项往往可以用非随机的函数进行刻画，剩下的随机噪声项，往往会具有某种平稳波动性，即在某条直线上下跳跃。平稳序列是用来描述某一种具有平稳波动性序列的序列，其定义如下。

如果时间序列$\{X_t\}$满足：

$\forall t\in \N,{\rm E}X_t^2<\infty$，即二阶矩存在；
$\forall \in\N,{\rm E}X_t=\mu$，即均值一致；
$\forall t,s\in\N,{\rm E}[(X_t-\mu)(X_s-\mu)]=\gamma_{t-s}$，即自协方差只与时间差有关。

就称$\{X_t\}$是平稳时间序列，称$\{\gamma_t\}$为$\{X_t\}$的自协方差函数。

从平稳序列的定义可以看出，它的平稳表现在两个方面，一是均值、方差的平稳性，即均值、方差与时间无关；二是相关性的平稳性，即序列中的任意两个随机变量自协方差函数，只与时间差有关，而与它们的绝对位置无关。
需要注意，自协方差函数是包含分布的方差的，因为${\rm D}X_t=\gamma_0$。这也说明了方差与时间无关，因为对任何$X_t$，其方差都是$\gamma_0$，是一个常数。如果$\gamma_0=0$，那么随机变量就是一个常数，没有讨论的必要，因此我们总假定$\gamma_0>0$。
从平稳序列的定义来看，它最重要的元素无疑是自协方差函数，这刻画了序列内部的关系。首先，很显然对于任何实时间序列，其自协方差序列都是实数列，除此外自协方差函数有以下三条重要性质：

1、对称性，即$\gamma_k=\gamma_{-k}$对所有$k\in\Z$成立。
2、非负定性，即对任何$n\in\N$，$n$阶自协方差矩阵

$\Gamma_n=\begin{bmatrix}
\gamma_0&\gamma_1&\cdots&\gamma_{n-1}\\
\gamma_1&\gamma_0&\cdots&\gamma_{n-2}\\
\vdots&\vdots&&\vdots\\
\gamma_{n-1}&\gamma_{n-2}&\cdots&\gamma_0
\end{bmatrix}$

总是非负定的。
3、有界性，即对任何$k\in\Z$，有$|\gamma_k|\le \gamma_0$。

同时满足以上三条性质的实数列称为非负定序列，平稳序列的自协方差函数就是非负定序列，并且可以证明，每个非负定序列都可以是一个平稳序列的自协方差函数。这里建立了非负定序列与平稳序列的对应性。
接下来对平稳序列的这三条性质进行证明。对称性最显然，由定义就可以直接看出，即

$\gamma_k={\rm Cov}(X_{t},X_{t+k})={\rm Cov}(X_{t+k},X_{t+k-k})=\gamma_{-k}.$

非负定性，即自协方差矩阵是非负定矩阵，对于任何一个二次型要证明其非负，就任取一个常数向量，计算

$\begin{aligned}
\boldsymbol a'_n\Gamma_n\boldsymbol a_n=&\sum_{j=1}^n\sum_{k=1}^na_ja_k\gamma_{j-k}\\
=&\sum_{j=1}^n\sum_{k=1}^na_ka_j{\rm Cov}(X_k,X_{j})\\
=&{\rm D}(\sum_{j=1}^n a_jX_j)\ge0.
\end{aligned}$

这个证明过程中，需要牢记的是二次型的写法，即将二次型写成一个双重求和的结果，每一项是$b_jb_k$与二次型矩阵的第$(j,k)$项乘积；并且将双重求和转化成一个单次求和的函数，这个思想也很重要。
有界性，用到柯西不等式，将随机变量中心化，即$Y_t=X_t-\mu$，那么${\rm D}Y_t={\rm D}X_t=\gamma_0$，${\rm Cov}(X_t,X_{t+k})={\rm Cov}(Y_t,Y_{t+k})=\gamma_k$，就有$|\gamma_k|=|{\rm E}(Y_tY_{t+k})|\le \sqrt{{\rm E}Y_t^2{\rm E}Y_{t+k}^2}=\gamma_0$，这里小于等于号就是柯西不等式的结果。
由有界性可以知道$-1\le \gamma_k/\gamma_0\le 1$，与相关系数有很大的相似之处，所以我们将$\gamma_k/\gamma_0$定义为平稳序列的自相关系数，也就是自协方差函数的归一化。在某些情况下，自相关函数甚至比自协方差函数还要重要。
我们再将目光投射到三条性质中，最不平凡的那条，即非负定性上。既然我们知道$\forall \boldsymbol a_n$，有$\boldsymbol a_n'\Gamma_n\boldsymbol a_n\ge 0$，那么作为临界情况的等号成立时意味着什么呢？显然等号很难对于所有$\boldsymbol a_n$都成立（除非$\Gamma_n=O$，但这是没有意义的），所以我们讨论对某个特定的$\boldsymbol a_n$等号成立的情况。由于

$\boldsymbol a_n'\Gamma_n\boldsymbol a_n={\rm D}(\sum_{j=1}^n a_jX_j),$

我们不妨定义$\boldsymbol X=(X_1,\cdots,X_n)'$，那么$\boldsymbol a_n'\Gamma_n\boldsymbol a_n={\rm D}(\boldsymbol a'\boldsymbol X)=0$，也就说明$\boldsymbol a'\boldsymbol X$是常数，结合其均值来看应该有$\boldsymbol a'X=\mu\boldsymbol a'\boldsymbol 1_n$（$\boldsymbol 1_n$指全是1的列向量）。由于我们规定$\boldsymbol a_n\ne0$，那么一定存在一个下标最大的分量$a_k\ne 0$，使得$X_k$可以被$X_1,\cdots,X_{k-1}$线性表示。这时，我们称$X_1,\cdots,X_n$是线性相关的。
并且进一步看，由于自协方差函数与序列位置无关，即

$\boldsymbol a_n'\Gamma_n\boldsymbol a_n=\sum_{j=1}^n\sum_{k=1}^na_ja_k{\rm Cov}(X_{t+j},X_{t+k})={\rm D}(\sum_{j=1}^n a_jX_{t+j})=0,$

所以对任何一组的连续的$(X_{t+1},\cdots,X_{t+n})$，都有$X_{t+k}$可以被$X_{t+1},\cdots,X_{t+k-1}$线性表示，并且表示系数是相同的。这一性质，表明对于退化的$\Gamma_n$，任何$X_{t},t\ge n$都可以被$X_0,\cdots,X_{n-1}$线性表示，这进一步说明了对于任意的$n$个$X_t$，它们一定是线性相关的，不管是不是连续增长的时间指标$t$。
事实上，用多元统计的观点看，设$\boldsymbol X=(X_{t+1},X_{t+2},\cdots,X_{t+n})$，那么$\Gamma_n={\rm D}(\boldsymbol X)$，即随机向量的协方差矩阵，那么自然有

${\rm E}(A\boldsymbol X+B)=A\boldsymbol X+B,\quad {\rm D}(A\boldsymbol X+B)=A{\rm D}(\boldsymbol X)A'=A\Gamma_nA'.$

当$A=\boldsymbol a'_n,B=c$的时候，显然有$D(\boldsymbol a_n'\boldsymbol X+c)=\boldsymbol a_n'\Gamma_n\boldsymbol a_n\ge 0$。
需要注意的是，平稳序列并不一定是平稳但散乱的，也可以具有很强的周期性，其典型例子就是调和平稳序列$X_t=b\cos(at+U),U\sim U(-\pi,\pi)$，它的自协方差函数是$\frac 12b^2\cos ((t-s)a)$，具有很强的周期性，所以观测样本也会具有周期性。
3.特殊的平稳序列——白噪声
白噪声是一种最为简单，但也颇具地位的平稳序列，其定义如下。

设$\{\varepsilon_t\}$是一个平稳序列，如果对任何$s,t\in\N$，都有

${\rm E}\varepsilon_t=\mu,\quad {\rm D}\varepsilon_t=\sigma^2,\\
\quad {\rm Cov}(\varepsilon_t,\varepsilon_s)=0,\quad t\ne s.$

就称$\{\varepsilon_t\}$是一个白噪声，记作${\rm WN}(\mu,\sigma^2)$。

关于其方差和协方差的另一种写法是

${\rm Cov}(\varepsilon_t,\varepsilon_s)=\left\{
\begin{array}l
\sigma^2,&t=s\\
0,&t\ne s
\end{array}=\sigma^2\delta_{t-s}.
\right.$

这里$\delta_k$是克罗内克(Kronecker)函数，当$k=0$时$\delta_k=1$，否则$\delta_k=0$。
白噪声又可以细分为以下几类：

当$\{\varepsilon_t\}$是独立序列时，称为独立白噪声（定义只保证了不相关）；
当$\mu=0$时，称为零均值白噪声；
当$\mu=0,\sigma^2=1$时，称为标准白噪声；
当$\{\varepsilon_t\}$服从正态分布且是独立序列时，称为正态白噪声。

4.多平稳序列的相互关系
多平稳序列的相互关系，指的是对于两个平稳序列$\{X_t\}$和$\{Y_t\}$，它们之间具有的相互性质。具体可以细分为正交平稳序列和不相关平稳序列，其定义如下：

正交的：如果$\forall s,t\in\Z$，都有${\rm E}(X_tY_s)=0$，就称$\{X_t\},\{Y_t\}$是正交的。
不相关的：如果$\forall s,t\in\Z$，都有${\rm E}(X_tY_s)={\rm E}X_t{\rm E}Y_s$，就称$\{X_t\},\{Y_s\}$是不相关的。

这两个定义很好从字面意义上理解。正交是垂直的推广，在线性代数中两个向量$a,b$正交被定义为其内积$\langle a,b\rangle=0$，在平稳序列中，就是乘积的期望为0；不相关就是二者不对对方产生影响，所以乘起来求期望与分开求期望相乘得到的结果理应是一样的。如果${\rm E}X_t{\rm E}Y_s=0$，那么正交序列和不相关序列本身等价，也就是说，对于零均值平稳序列，其正交性和不相关性是等价的。
为什么要讨论这两种特殊的关系呢？我们以后可能会对平稳序列进行求和，即$Z_t=X_t+Y_t$，如果$\{Z_t\}$本身也能够是平稳序列那再好不过了。幸运的是，$\{X_t\},\{Y_t\}$是正交、不相关序列时，都能让$\{Z_t\}$是平稳序列。
要证明$\{Z_t\}$是平稳序列，就要证明其二阶矩有限、期望平稳、自协方差函数仅与时间差有关。期望平稳是显然的，有$\mu_Z=\mu_X+\mu_Y$；二阶矩有限也是显然的，有

${\rm E}Z_t^2={\rm E}(X_t+Y_t)^2\le 2{\rm E}X_t^2+2{\rm E}Y_t^2<\infty.$

接下来对正交、不相关序列，分别求$\{Z_t\}$的自协方差函数。首先是正交的情况，有

$\begin{aligned}
{\rm Cov}(Z_{t},Z_s)=&{\rm Cov}(X_t,X_s)+{\rm Cov}(X_t,Y_s)+{\rm Cov}(Y_t,X_s)+{\rm Cov}(Y_t,Y_s)\\
=&\gamma_X(t-s)+\gamma_Y(t-s)+{\rm E}(X_sY_t)-{\rm E}X_s{\rm E}Y_t+{\rm E}(X_tY_s)-{\rm E}X_t{\rm E}Y_s\\
=&\gamma_X(t-s)+\gamma_Y(t-s)-2\mu_X\mu_Y;
\end{aligned}$

然后是不相关的情况，立马得到${\rm Cov}(Z_t,Z_s)=\gamma_X(t-s)+\gamma_Y(t-s)$。这两个数都是$t-s$的函数，这就证明了对正交、不相关平稳序列，其和仍然是平稳序列。
加和的自协方差函数不方便记忆，可以记以下的简化结论：对于零均值的正交平稳序列$\{X_t\},\{Y_t\}$，他们的和$\{Z_t\},Z_t=X_t+Y_t$仍是平稳序列，且$\mu_Z=\mu_X+\mu_Y,\gamma_Z(k)=\gamma_X(k)+\gamma_Y(k)$。
回顾总结

任何时间序列经过适当的变换，都可以拆解为趋势项、季节项、随机噪声，并且趋势项和季节项一般被认为是非随机函数。
平稳序列是二阶矩存在、期望一致、自协方差只与时间差有关的时间序列，满足这三个条件就是平稳序列，这一般被用来验证序列的平稳性。
平稳序列中最重要的是自协方差函数$\{\gamma_k\}$，这是一个实数列，满足对称性、非负定性、有界性三个性质。
满足对称性、非负定性、有界性的实数列被称为非负定序列，一个非负定序列一定是某个平稳序列的自协方差函数。但非负定性的验证比较麻烦，所以验证一个序列是非负定序列一般是证明它是某个平稳序列的自协方差函数。
如果$\Gamma_n$退化，则$X_1,\cdots,X_n$线性相关，并且可以证明任何下标不小于$n$的项$X_t$都可以用$X_1$到$X_{n-1}$这$n-1$项线性表示。
平稳序列是特殊的时间序列，白噪声${\rm WN}(\mu,\sigma^2)$是特殊的平稳序列，它的主要特征是序列不相关性，也就是序列之间任意两个不同的随机变量无关，一样需要满足均值、方差的一致性。
白噪声中，又有独立白噪声、零均值白噪声、标准白噪声、正态白噪声几类特殊白噪声。
平稳序列正交指${\rm E}(X_tY_s)=0$，不相关指${\rm E}(X_tY_s)={\rm E}X_t{\rm E}Y_s$，对零均值平稳序列这两个定义是等价的。
平稳的正交、不相关序列加和仍然是平稳序列，且对于零均值的情况，自协方差函数为两个分开的自协方差函数之和。

有些时间序列具有非常显著的趋势，有时我们分析的目的就是要找到序列中的这种趋势，并利用这种趋势对未来的发展作出合理的预测。



趋势拟合法

趋势拟合法就是把时间作为自变量，相应的序列观测值做为因变量，建立序列值随时间变化的回归模型的方法。根据序列所表现出的线性或者非线性特征，拟合方法又可以具体分为线性拟合曲线拟合。



线性拟合

如果长期趋势呈现出线性特征，那我们可以用线性模型来拟合。 

R语言使用lm函数拟合线性趋势。 

lm(Y ~ a + X1 + X2 +…+Xn,data = )


  
-Y:响应变量 

  -a:指定是否需要常数项 

  （1）a = 1，模型有非零常数项，这是默认设置。 

  （2）a = 0，模型不需要常数项。 

  -X1,…,Xn:自变量。 

  -data:数据框名。如果自变量和响应变量不是独立输入变量而是共同存在于某个数据框中，则需要指定数据框名。




#读入数据
x<-c(8444,9215,8879,8990,8115,9457,8590,9294,8997,9574,9051,9724,9120,
     + 10143,9746,10074,9578,10817,10116,10779,9901,11266,10686,10961,10121,
     + 11333,10677,11325,10698,11624,11502,11393,10609,12077,11376,11777,
     + 11225,12231,11884,12109)
#构造时间变量
t<-c(1:40)
#拟合回归模型
x.fit<-lm(x~t)
#查看拟合信息
summary(x.fit)
Call:
lm(formula = x ~ t)

Residuals:
    Min      1Q  Median      3Q     Max 
-853.06 -329.64   63.54  314.46  794.04 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8491.765    137.964   61.55   <2e-16 ***
t             90.009      5.864   15.35   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 428.1 on 38 degrees of freedom
Multiple R-squared:  0.8611,    Adjusted R-squared:  0.8575 
F-statistic: 235.6 on 1 and 38 DF,  p-value: < 2.2e-16




#绘制拟合效果图
x<-ts(x)
plot(x)





abline(lm(x~t),col=2)





曲线拟合

如果长期趋势呈现出非线性特征，那么我们可以用曲线典线模型来拟合它。对曲线模型进行参数估计时，能转换成线性模型的都转化成线性模型，用线性最小二乘法进行参数估计，实在不能转化成线性模型的，就用迭代法进行参数估计。 



R语言针对非线性趋势的拟合也分为两类：一类可以写成关于时间t的多项式，这时仍然可以用lm函数拟合。另一类无法通过适当的变换变成线性回归模型，只能通过非线性回归解决，这时要用nls函数。 

nls(Y ~ f(x1,…,xn),data=,start=)


  
-Y:响应变量。 

  -X1,…,Xn：自变量 

  -f:非线性函数 

  -data:数据框名 

  -start:如果需要利用迭代计算未知参数，可以指定迭代初始值。


1949-2008年化肥产量序列进行曲线拟合’ 

lm函数拟合



a<-read.table("D:/R-TT/book4/4R/data/file12.csv",sep=",",header = T)
x<-ts(a$output,start=1949)
#lm函数拟合
t1<-c(1:60)
t2<-t1^2
x.fit1<-lm(x~t1+t2)
summary(x.fit1)
Call:
lm(formula = x ~ t1 + t2)

Residuals:
    Min      1Q  Median      3Q     Max 
-532.12 -164.92   24.68  105.51  716.37 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 319.0255   105.3371   3.029  0.00369 ** 
t1          -57.7690     7.9679  -7.250 1.22e-09 ***
t2            2.3551     0.1266  18.601  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 263 on 57 degrees of freedom
Multiple R-squared:  0.9755,    Adjusted R-squared:  0.9746 
F-statistic:  1133 on 2 and 57 DF,  p-value: < 2.2e-16

nls函数拟合



 x.fit2<-nls(x~a+b*t1+c*t1^2,start = list(a=1,b=1,c=1))
> summary(x.fit2)

Formula: x ~ a + b * t1 + c * t1^2

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 319.0255   105.3371   3.029  0.00369 ** 
b -57.7690     7.9679  -7.250 1.22e-09 ***
c   2.3551     0.1266  18.601  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 263 on 57 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.059e-08

根据输出结果可以知道，这两个函数得到的拟合结果完全一致，1949-2008年我国化肥序列的拟合模型为： 

Xt = 319.0255 - 57.769t +2.3551t^2     N（0，263^2）



y<-predict(x.fit2)
y<-ts(y,start = 1949)
plot(x,type = "p")
lines(y,col=2,lwd=2)





平滑法

平滑法是进行趋势分析和预测时常用的一种方法，它是利用修匀技术，消弱短期随机波动的序列的影响，使序列平滑化，从而显示出变化的规律。它具有调节灵活、计算简便的特征，广泛应用于计量经济、人口研究等诸多领域。根据所有的平滑技术的不同，平滑法又可以具体分为移动平均法与指数平滑法。



移动平均法

移动平均法的基本思想是对于一个时间序列{Xt},我们可以假定在一个比较短的时间间隔内，序列的预测值是比较平稳的，它们之间的差异主要是由随机波动造成的，根据这种假定，我们可以使用一定时间间隔的平均值作为下一期的估计值。 

移动平均的期数对原序列的修匀效果影响很大，要确定移动平均的期数，一般会从如下三个方面加以考虑： 

（1）事件的发展有无周期性。 

（2）对趋势平滑性的要求。 

（3）对趋势反映近期灵敏度的要求。 

在R语言中，TTR程序包中的SMA函数是专门用来作简单移动平均趋势拟合的函数。 

SMA(x,n)


  
-x:需要做简单的移动平均的序列名 

  -n:移动平均期数


对北京市1949年-1998年每年最高气温序列进行进行5期移动平均拟合。



library(TTR)
a<-read.table("D:/R-TT/book4/4R/data/file6.csv",",",header = T)
x<-ts(a$temp,start = 1949)
x.ma<-SMA(x,n=5)
plot(x,type = "o")
lines(x.ma,col=2,lwd=2)





指数平滑法

移动平均法实际上就是用一个简单的加权平均数作为某一期的趋势的估计值。 

原理为近期的结果对现在的影响大些，远期的结果对现在的影响小些。 

各期权重随时间间隔的增大而呈指数衰减，这就是指数平滑法的基本思想。 

1.简单的指数平滑 


2.Holt两参数指数平滑 

适用于对含有线性趋势的序列进行修匀。它的基本思想是假定序列有一个比较固定的线性趋势——每期都递增r或者递减r,那么第t期的估计值就应该等于第t-1期的观察值加上每期固定的趋势变动。 



3.Holt-Winters三参数指数平滑 

假定要进行指数平滑的序列为{Xt},{Xt}序列既含有趋势又含有季节。这个季节因子可以随每年的具体情况波动，因为引入了季节因子，所以构建了三参数模型。 

在R语言中，HoltWinters函数可以完成上述三种平滑趋势拟合。 

HoltWinters(x,alpha=,beta=,gamma=,seasonal=)


  
-x:要进行指数平滑的序列名 

  -alpha:随机波动部分的参数 

  -beta:趋势部分的参数 

  -gamma:季节部分的参数 

  这三个指数联合起来，确定要拟合的指数平滑模型类型： 

  （1）当alpha不指定时，beta = F,gamma=F,表示拟合简单指数平滑模型。 

  （2）当alpha和beta不指定时，gamma=F时，表示拟合Holt两参数指数平滑模型。 

  （3）当三个参数都不指定时，表示拟合Holt-Winters三参数指数平滑模型。 

  -seasonal:当既含有季节的又含有趋势时，指定季节与趋势的关系。 

  seasonal = “additive”表示加法关系， 这是系统默认选项。seasonal=”multiplicative”表示乘法关系。


对1964-1999年中国纱年产量序列分析,并预测未来10年的



#读入序列
a<-read.table("D:/R-TT/book4/4R/data/file4.csv",sep=",",header = T)
x<-ts(a$output,start = 1964)
#进行Holt两参数平滑
x.fit<-HoltWinters(x,gamma = F)
x.fit
Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
HoltWinters(x = x, gamma = F)

Smoothing parameters:
 alpha: 0.855644
 beta : 0.158537
 gamma: FALSE

Coefficients:
       [,1]
a 565.55301
b  12.29066



#绘制Holt两参数指数平滑拟合效果图
plot(x.fit)





#预测序列并绘制预测效果图
library(forecast)
x.fore<-forecast(x.fit,h=10)
x.fore
 Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2000       577.8437 545.8643 609.8230 528.9355 626.7519
2001       590.1343 545.1051 635.1636 521.2681 659.0006
2002       602.4250 544.7497 660.1003 514.2182 690.6318
2003       614.7157 544.3115 685.1198 507.0417 722.3896
2004       627.0063 543.6025 710.4101 499.4512 754.5614
2005       639.2970 542.5375 736.0565 491.3161 787.2779
2006       651.5876 541.0751 762.1002 482.5733 820.6020
2007       663.8783 539.1960 788.5606 473.1931 854.5635
2008       676.1690 536.8922 815.4458 463.1635 889.1744
2009       688.4596 534.1622 842.7570 452.4821 924.4371
plot(x.fore)



对1962年1月和1975年12月平均每头奶牛月产奶量序列进行Holt-Winters三参数指数平滑



#读入序列
b<-read.table("D:/R-TT/book4/4R/data/file5.csv",sep=",",header = T)
x<-ts(b$milk,start = c(1962,1),frequency = 12)
#进行Holt-winters三参数指数平滑
x.fit<-HoltWinters(x)
x.fit
Holt-Winters exponential smoothing with trend and additive seasonal component.

Call:
HoltWinters(x = x)

Smoothing parameters:
 alpha: 0.68933
 beta : 0
 gamma: 0.8362592

Coefficients:
          [,1]
a   885.775547
b     1.278118
s1  -16.743296
s2  -59.730034
s3   47.492731
s4   56.203890
s5  115.537545
s6   84.554817
s7   39.580306
s8   -4.702033
s9  -54.554684
s10 -51.582594
s11 -85.953466
s12 -42.907363

#绘制Holt-winters三参数指数平滑拟合效果图
plot(x.fit)





#预测序列并绘制预测效果图
x.fore<-forecast(x.fit,h=24)
plot(x.fore)




月度产奶量Holt-Winters三参数指数平滑序列预测图