一、基础知识

时间序列分析就是对一个时间序列进行建模，扣除各种趋势项（线性趋势、余弦趋势、有色噪声ARIMA），得到一个白噪声序列；换言之，我们要提取其中的有用信息（非白噪声序列），最后剩下的越接近白噪声越说明提取的越号，越充分。

基础知识部分整理了几个子博客，链接如下：
【时序分析】平稳性、可逆性与几个重要过程
【时序分析】自相关函数与偏自相关函数（R语言）
【时序分析】滑动平均过程 MA
【时序分析】自回归过程 AR
【时序分析】自回归滑动平均混合模型 ARMA

二、例题

1、3.5（P37）

数据文件 wages包含了1981年7月到1987年6月美国服装和纺织品行业工人的平均时薪（以美元计)的月度值.

• (a)画出并解释这些数据的时间序列图.
• (b)对该时间序列用最小二乘法拟合线性时间趋势．解释回归结果．保存拟合模型的标准残差以便进一步分析.
• (c) 构造并解释来自(b)的标准残差的时间序列图.
• (d)对工资时间序列用最小二乘法估计二次时间趋势．解释回归结果．保存标准残差以便进一步分析.
• (e)绘出并解释(d)中标准残差的时间序列图.

【解】
第(a)问中需要先导入TSA数据库，否则无法data(wages)

library(TSA);
data(wages);
plot(wages,ylab='Monthly Wages',type='o');

第(b)(c)问

wages.lm=lm(wages~time(wages));
summary(wages,lm);
y=rstudent(wages.lm);
plot(y,x=as.vector(time(wages)),ylab='Standardized Residuals',type='o');

第(d)(e)问

wages.lm2=lm(wages~time(wages)+I(time(wages)^2));
summary(wages.lm2);
y=rstudent(wages.lm2);
plot(y,x=as.vector(time(wages)),ylab='Standardized Residuals',type='o');

2、自相关和偏相关

计算并画出以下AR(2)模型的自相关函数 ${\varphi }_1=-0.4,{\varphi }_2=0.5$：

Y=ARMAacf(ar=c(-0.4,0.5),lag.max=12)
Y
acf(Y)
___
1      2      3      4      5      6      7      8      9     10     11     12     13
-0.238 -0.142 -0.002 -0.002 -0.003 -0.003 -0.004 -0.004 -0.005 -0.005 -0.006 -0.006 -0.007

3、第一次综合作业

请以卫星重力梯度观测值 Vzz 分量（重力位二阶径向导数） 的时间序列（见数据文件“data.txt”）为例，进行如下分析：
（1） 画出并解释这些数据的时间序列图。判断是否具有相关性？序列是否具有趋势项或周期项？
（2） 如果具有趋势项，请对该时间序列用最小二乘法拟合线性时间趋势。并画出扣除线性趋势后的时间序列图。
（3） 分析（2）中扣除线性趋势后的时间序列图，如果具有周期性，请用余弦趋势模型进行拟合， 并画出扣除余弦趋势后的时间序列图。
（4） 构造并解释来自（3）的标准残差（学生化残差） 的时间序列图，分析其是否满足正态性？（QQ 图，直方图）
（5） 如果并不满足正态性，是否可采用 ARMA（p，q）模型对来自（3） 中的时间序列（即扣除了线性趋势（2） 和余弦趋势（3） 的残差时间序列） 进行建模？ 请绘制自相关函数和偏自相关函数图， 能否对 ARMA（p，q）模型定阶？
数据说明： data.txt 文件中包含 2 列，第 1 列为时刻（单位： s），第2 列为 Vzz 梯度值（单位： mE） ，观测值间隔为 30s。

第一、二问：

library(TSA);
f1='C:/gwork/study/postgraduate/时间序列/练习/data.txt'
data<-read.table(f1)
Vzz<-ts(data[2])
time<-ts(data[1])
model=lm(Vzz~time)					# 最小二乘拟合
plot(time,Vzz,type='o',
ylab=expression(Vzz),
xlab=expression(SOD/s),col='red')	# 绘制原始数据
abline(model,col='blue')			# 绘制拟合线
Vzz1<-ts(residuals(model))			# 取出拟合残差
plot(time,Vzz1)						# 绘制拟合残差

第三、四问：

rs<-ts(as.vector(Vzz1),freq=180)        # 将残差1构造成时间序列
har.=harmonic(rs,1)    
model2=lm(rs~har.)                      # 余弦趋势拟合
Vzz2<-ts(residuals(model2));            # 将残差取出
plot(ts(fitted(model2),freq=180),x=ts(time,freq=180),type='l',
ylab=expression(Vzz1),
xlab=expression(SOD/s),col='blue')      # 绘制余弦拟合结果
plot(time,Vzz2)                         # 绘制余弦拟合残差
qqnorm(Vzz2)                            # 绘制QQ图，若为直线则为正态，否则不是

第五问：

acf(Vzz2)                               # 自相关函数
pacf(Vzz2)                              # 偏相关函数
library(zoo)
library(forecast)
auto.arima(Vzz2)