时间序列分析
 
需要解释的概念
 
平稳时间序列
 
平稳时间序列
 
平稳时间序列及其数字特征
 
 
均值为常数，协方差函数与时间无关的时间序列叫做平稳时间序列
 
 
平稳时间序列有一种简化模型，就是均值为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。
 
 
大概总结一下，如下图所示
 
 

童鞋们觉得文章不错，就麻烦点一下下面人工智能的教程链接吧，然后随便翻阅一下
 
https://www.captainbed.net/qtlyx
 
 
 
1.非平稳时间序列
 
之前我们说明了怎么样的时间序列是序列平稳的，但是世界并不是那么美好，很多时间序列都不是平稳序列，所以这里就要求我们做一些处理了。
 
首先我们来看一下非平稳时间序列长什么样。在AR模型中，只要自回归系数都绝对值都是小于1的，那么序列就是平稳的，所以这样一个序列，自回归系数等于1，就是不平稳的序列了。
 
yt = yt-1 + c
 
c是一个服从正态分布的噪音。
 
 
 
#example 10
set.seed(12345)
ut = rnorm(50,0,1.5)
xt = cumsum(ut)
plot(xt,type = 'o');abline(h = 0)
 
其中，cumsum是一个计算累计数的函数。比如cumsum(c(1,2,3,4,5))=(1,1+2,1+2+3,1+2+3+4.....)=(1,3,6,10......)
 
 
 

  
 
 
 
这就是对非平稳序列的一个直观的感受了。
 
2.非平稳序列的平稳方法--差分
 
非平稳序列往往一次到两次差分之后，就会变成平稳序列。什么是差分呢？就是后一时间点的值减去当前时间点，也就是yt-yt-1。
 
 
 
#example 11
x = 1:10
diff(x,d=1)
diff(x,d=2)
 
 
 
这里，我们对1,2,3,4,5......这个序列做了两次差分，都是后一个数减去前一个数。
 
 
值得注意的一点是，每一次差分之后，都会少一个序列值。
 
 
 
#example 12
plot(diff(xt,d = 1),type = 'o');abline(h = 0)
plot(diff(xt,d = 2),type = 'o');abline(h = 0)
 
我们用之前的序列试一下水，可以看到，一阶差分和二阶差分后，看上去都平稳了呢！
 
 
 

  
 
 
 
 
 
 
 
 
 
 
3.判断序列是否平稳
 
 
 
前面我们用肉眼看了序列的平稳性，但是作为一个时间序列分析者，竟然用眼睛主观判断，这有点不合逻辑。很幸运的是，我们根据单位根过程（有兴趣的读者查找相关资料），可以进行假设检验，譬如DF与ADF检验。
 
adf检验函数包含在tseries这个包中，使用前我们要先引用它。
 
 
 
#example 13
adf.test(xt)
adf.test(diff(xt))
adf.test(diff(xt,d=2))
 
 
 
 
大家注意看哦，当没有做差分的时候，p-value是0.47+，而备择假设是stationary，也就是平稳时间序列，所以零假设就是非平稳时间序列。p-value>0.05的时候，在95%的置信度下，我们是不能拒绝原假设的，所以我们不能说xt原序列是时序平稳的，但是对于一阶差分和二阶差分就是可以的了。
 
 
 
童鞋们觉得文章不错，就麻烦点一下下面人工智能的教程链接吧，然后随便翻阅一下
 
https://www.captainbed.net/qtlyx

 
时间序列分析：平稳时间序列分析
 
一.录入数据
1.读入数据
2.转为时间序列数据
  
二. 平稳时间序列建模
1.建模步骤
2.根据序列的时序图，自相关图和偏自相关图判断序列的平稳性，随机性
2.1.时序图检验
2.2.自相关图检验
2.3判断该序列的纯随机性
   
3.选择适当模型
4.模型检验及优化
4.1模型显著性检验
4.2参数显著性检验
4.3模型优化
  
  
三.时间序列预测

 
一.录入数据
 
1.读入数据
 
在R中输入一下命令:
 
rain=read.table("/home/ddy/桌面/习题3.17数据.txt",sep="\t",header = F)
> head(rain)
 

 我们发现该数据录入的时候是按行录入的，我们需要的时按列录入的数据，所以我们对原数据先转置，再将数据的每一列链接到第一列下面，代码实现如下：
 
#将数据转置
> rain=t(rain)
> rain_m=matrix(data=0, nrow = 64, ncol = 1)
> k=1
#将数据的每一列链接到第一列下面
> for(i in 1:8){  
   for(j in 1:8){ 
      rain_m[k]=rain[j,i]  
      k=k+1   
     }
}
> rain_m
 
2.转为时间序列数据
 
rain_m=ts(rain_m)
 
二. 平稳时间序列建模
 
1.建模步骤
 
 
2.根据序列的时序图，自相关图和偏自相关图判断序列的平稳性，随机性
 

 
 
 
2.1.时序图检验
 
由时序图可初步得知该地区年降雨量几乎都在80在右波动，没有明显的周期性和单调性，所以认为该序列为平稳序列。
 
2.2.自相关图检验
 
在自相关图中，我们发现除了延迟1,2阶的自相关系数在2倍标准差以外，其他阶数的自相关系数2倍标准差以内波动，根据自相关系数的这个特点可以判断该序列具有短期相关性，即认为该序列为平稳序列。
 
并且非0自相关系数衰减到小值是一个突然的过程，不是一个连续渐变的过程，这是自相关系数截尾的的典型特征。
 
在偏自相关系数衰减到0的过程比较缓慢，可知道偏自相关系数不截尾。
 
综合1,2我们可认为该地区年降雨量序列是平稳序列。
 
2.3判断该序列的纯随机性
 
 
根据检验结果我们可以得出，从延迟6阶的检验结果中我们可以把得到P值都小于显著性水平0.05。因此，我们拒绝序列纯随机的原假设，即该序列为非白噪声序列，我们可以认为该该地区年降雪量序列变动不属于随机变动，这说明我们可以根据历史信息来预测未来年份的该地区年降雪量。
 
3.选择适当模型
 
在自相关系数衰减到0的过程，除了延迟2阶的偏自相关系数在2倍标准差以外，其他阶数的偏自相关系数2倍标准差以内波动，这是一个自相关系数2阶拖尾的的典型特征。
 
所以我们选择初步用MA(2) 模型来拟合序列，为了尽量避免因个人经验不足而导致模型识别不准确，所以我们运用auto.arima()函数来自动识别模型阶数，并从中得到参数估计值。
 auto.arima(rain_m)的结果为：
 
 由此可知，我们不能选择MA(2)模型，而要选择MA(1)模型，接下来我们对MA(1)模型进行检验：
 
 
4.模型检验及优化
 
4.1模型显著性检验
 
模型显著性检验只要是检验模型的有效性，一个模型是否显著有效主要看它提取的信息是否充分，一个好的拟合模型应该能够提取观察值序列中几乎所有的样本相关信息。换一句话说就是拟合残差项中将不在蕴含任何相关信息，即残差序列应该为白噪声序列，这样的模型才称为显著有效模型。所以我们可以认为模型的显著性检验就是残差序列的白噪声简检验。
 

 由检验结果可知，各阶延迟下的LB统计量的P值都大于0.05,可以认为这个拟合模型的残差序列属于白噪声序列，即该拟合模型显著有效。
 
此外，即使通过了MA(1)模型的检验。为了充分考虑各种可能，我们建立多个拟合模型，我们发现MA(2)模型的结果并没有 MA(1)的结果显著，所以我们最终选择MA(1)模型进行拟合。
 
4.2参数显著性检验
 
检验函数：
pt(t,df = ,lower.tail = )

t:t统计量的值图（该例子中：t=0.2103/0.0982）
df:自由度（n-参数个数）
lower.tail:确定计算概率的方向。对于参数显著性检验，如果参数估计值为负，选择lower.tail=T
 
4.3模型优化
 
若一个模型通过里检验，说明在一定置信水平下，该模型能够有效的拟合观察值序列的波动，但这种模型并不是唯一的。此时我们可以通过比较各个模型的AIC值和SBC值来进行最优模型的选择。
 
三.时间序列预测
 
利用拟合模型，预测城市未来5年的降雪量：
 
> rain.fore=forecast(rain_m_fit,h=5)
 
 
> plot(rain.fore)
 
得到该城市未来5年的降雪量预测图：
 
 
代码汇总：
rain=read.table("/home/ddy/桌面/习题3.17数据.txt",sep="\t",header = F)
head(rain)

rain=t(rain)

rain_m=matrix(data=0, nrow = 64, ncol = 1)
k=1
for(i in 1:8){
  for(j in 1:8){
    rain_m[k]=rain[j,i]
    k=k+1
  }
}

rain_m=data.frame(rain_m[-64,])
head(rain_m)

rain_m=ts(rain_m)
par(mfcol=c(3,1))
plot(rain_m,main = "年降雪量时序图",col=6)
acf(rain_m,main = "年降雪量自相关图",col=2)
pacf(rain_m,main = "年降雪量偏自相关图",col=2)

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

arima(x = rain_m, order = c(0, 0, 1), method = "ML")
for(i in 1:2)print(Box.test(rain_m_fit$residuals,lag=6*i))
pt( 0.2103/0.0982,df=62,lower.tail=F)
rain.fore=forecast(rain_m_fit,h=5)
rain.fore
plot(rain.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(L2,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(U2,col=4,lty=2)

 
文章目录
 
前言
一、方差不平稳
二、具有趋势性
1、确定性趋势的删除——趋势拟合法
2、随机趋势的删除——差分法
3、未知趋势的删除——平滑法
1. 移动平均法
2. 指数平滑法
   
4、季节性
1.Method1: 减掉趋势成分后取平均
2.Method2：使用哑元变量做回归分析
  
 

 
 
前言
 
时间序列为什么要平稳后建模？
 
时间序列是平稳的，则没有任何预测价值，因为均值不变，所以实际上人们只会对非平稳时间序列感兴趣
之所以研究平稳时间序列，因为总希望可以通过差分等方式把非平稳序列变成平稳的，而平稳时间序列可以通过数学公式予以解释
 
 
 
一、方差不平稳
 
最简单的是做对数变换，使得前后方差没有那么变化。
 
#Log transform: variance stationary
data(AirPassengers)
AirPassengers
par(mfrow=c(1,2))
plot(AirPassengers)
plot(log(AirPassengers))
 
 
二、具有趋势性
 
1、确定性趋势的删除——趋势拟合法
 
最简单的处理方法就是考虑均值函数可以由一个时间的确定性函数来描述，比如可以用回归模型来描述。对于一个有线性趋势的序列，可以做一个线性模型。
 
Linear trend
 
##read the gasline data
##Abraham and Ledolter (1983) on the
##monthly gasoline demand in Ontario over the period 1960 - 1975.
gas = scan("gas.dat")
## note that read.table() would give an error here since the data is not in rectangular format.
gas.ts = ts(gas, frequency = 12, start = 1960)
gas.ts
plot(gas.ts, main = "Gasoline demand in Ontario", ylab = "Million gallons")
#fit trend:OLS with constant and trend
time = 1:length(gas.ts)
fit = lm(gas.ts~time)
summary(fit)#R-squared:  0.8225
par(mfrow=c(2,1))
plot(time, gas.ts, type='l', xlab='months', ylab='Million gallons gasoline',
main="OLS fit")
legend(.1, 256000, c("Time plot", "OLS fit"), col = c(1,2), text.col = "black",
merge = TRUE, lty=c(1,1))
abline(fit,col=2,lwd=2)
plot(time, fit$resid, type='l', xlab="months", ylab="Residuals",main="Residuals OLS [yt-fit(yt)]")
abline(a=0, b=0)#残差在0附近点波动
#savePlot("plot of Residuals_gas",type="pdf")
 
 
Quadratic trend
 
fit2 = lm(gas.ts~time+I(time^2))
summary(fit2)#R-squared:  0.8327
names(fit2)
par(mfrow=c(2,1))
plot(time, gas.ts, type='l', xlab='months', ylab='Million gallons gasoline',
     main="OLS fit")
legend(.1, 256000, c("Time plot", "OLS fit"), col = c(1,2), text.col = "black",
       merge = TRUE, lty=c(1,1))
lines(time,fitted(fit2),col=2)   
plot(time, fit2$resid, type='l', xlab="months", ylab="Residuals",main="Residuals OLS [yt-fit(yt)]")
abline(a=0, b=0)
 
 
Polynomial trend
 
2、随机趋势的删除——差分法
 
 随机趋势的删除——差分法。
 
一阶差分可以使线性趋势的序列实现趋势平稳
 
Zt = as.ts(rnorm(1000, sd = 20))
RW1 = as.ts(cumsum(Zt))
par(mfrow = c(2,2))
plot(RW1, main = "Random walk")
acf(RW1, main = "Correlogram random walk")
Yt=as.ts(diff(RW1))
plot(Yt,main='Diff(RW1)')
acf(Yt)
 
 
序列蕴含着曲线趋势，通常低阶（二阶或三阶）差分可提取出曲线趋势的影响 ，这里效果不好的话，可以先取对数，再做差分
对于蕴含着固定周期的序列进行步长为周期长度的差分运算，通常可以较好地提取周期信息
差分的阶数要适度，避免过差分
 
3、未知趋势的删除——平滑法
 
 平滑法——顾名思义，使序列变平滑。 期望结果，是显示出趋势变化的规律。
 
1. 移动平均法
 
假定在短的时间间隔内，序列的取值是比较稳定的，序列的大小差异主要是由随机波动造成的。用一段时间间隔内的平均值作为某一期的估计值。
 
简单移动平均
 
##Moving average Model
Vt5 = filter(Zt, rep(1/5, 5), sides = 2)
Vt21 = filter(Zt, rep(1/21, 21), sides = 2)
par(mfrow = c(2,2))
plot(Vt5, main = "Moving 5-point average of above white noise series")
acf(Vt5, main = "Correlogram 5-point moving average", lag.max = 40, na.action = na.pass)
plot(Vt21, main = "Moving 21-point average of above white noise series")
acf(Vt21, main = "Correlogram 21-point moving average", lag.max = 40, na.action = na.pass)
par(mfrow = c(1,1))
 
 
权重移动平均
 
2. 指数平滑法
 
4、季节性
 
Seasonality means an effect that happens at the same time and with the same magnitude and direction every year. 从而使得经过季节调整的序列能够较好的反应社会经济指标运行基本态势。
 
1.Method1: 减掉趋势成分后取平均
 
## Investigating seasonality
library(fields)
z = matrix(fit$resid, ncol=12, byrow=TRUE)
colnames(z) = c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
bplot(z, xlab="Month", ylab="Detrended demand of gasoline", main="Annual seasonality of the detrended demand of gasoline")    
mean=matrix(0, nrow=ncol(z), ncol=1) #mean for each month
for(i in 1:ncol(z)){
  mean[i,1]=mean(z[,i])
}
Seas = as.matrix(rep(mean,16)) # matrix with seasonal term
des.gas = fit$resid - Seas   # detrended and now deseasonized series
par(mfrow=c(2,1))
plot(fit$resid, type='l', xlab="months", ylab="Residuals", main="Residuals OLS [yt-fit(yt)]")
abline(a=0, b=0)
plot(des.gas, type='l', xlab="months", ylab="Deseasonized residuals", main="Deseasonized ResidualsOLS")
abline(a=0, b=0)
 
 
2.Method2：使用哑元变量做回归分析
 
##seasonal modeling and estimating seasonal trends
##such as for the average monthly temperature data 
library(TSA)  #install package TSA
win.graph(width=4.875, height=2.5,pointsize=8)
data(tempdub); 
plot(tempdub,ylab='Temperature',type='o')
month.=season(tempdub) # period added to improve table display
model2=lm(tempdub~month.) # -1 removes the intercept term
summary(model2)
 
Call:
lm(formula = tempdub ~ month.)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2750 -2.2479  0.1125  1.8896  9.8250 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       16.608      0.987  16.828  < 2e-16 ***
month.February     4.042      1.396   2.896  0.00443 ** 
month.March       15.867      1.396  11.368  < 2e-16 ***
month.April       29.917      1.396  21.434  < 2e-16 ***
month.May         41.483      1.396  29.721  < 2e-16 ***
month.June        50.892      1.396  36.461  < 2e-16 ***
month.July        55.108      1.396  39.482  < 2e-16 ***
month.August      52.725      1.396  37.775  < 2e-16 ***
month.September   44.417      1.396  31.822  < 2e-16 ***
month.October     34.367      1.396  24.622  < 2e-16 ***
month.November    20.042      1.396  14.359  < 2e-16 ***
month.December     7.033      1.396   5.039 1.51e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.419 on 132 degrees of freedom
Multiple R-squared:  0.9712,	Adjusted R-squared:  0.9688 
F-statistic: 405.1 on 11 and 132 DF,  p-value: < 2.2e-16